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fSJ , We study eccentric equatorial orbits of a test-body around a Kerr black hole under the influence 

' of gravitational radiation reaction. We have adopted a well established two-step approach: assuming 

' that the particle is moving along a geodesic (justiflable as long as the orbital evolution is adiabatic) 

we calculate numerically the fluxes of energy and angular momentum radiated to infinity and to 
?— j ' the black hole horizon, via the Teukolsky-Sasaki-Nakamura formalism. We can then infer the rate 

of change of orbital energy and angular momentum and thus the evolution of the orbit. The orbits 
are fully described by a semi-latus rectum p and an eccentricity e. We find that while, during the 
inspiral, e decreases until shortly before the orbit reaches the separatrix of stable bound orbits (which 
is deflned by ps(e)), in many astrophysically relevant cases the eccentricity will still be signiflcant in 
the last stages of the inspiral. In addition, when a critical value Pcrit(e) is reached, the eccentricity 
begins to increase as a result of continued radiation induced inspiral. The two values ps, Pcrit (for 
^ . given e) move closer to each other, in coordinate terms, as the black hole spin is increased, as they 

\^ ' do also for flxed spin and increasing eccentricity. Of particular interest are moderate and high 

QQ ^ eccentricity orbits around rapidly spinning black holes, with p(e) « Ps(e). We call these "zoom- 

. whirl" orbits, because of their characteristic behaviour involving several revolutions around the 

CO ' central body near periastron. Gravitational waveforms produced by such orbits are calculated and 

, shown to have a very particular signature. Such signals may well prove of considerable astrophysical 

importance for the future LISA detector. 

O 

O^' I. INTRODUCTION 

A. Background 

k>( \ Binary star systems, consisting of compact objects such as black holes and neutron stars, are relatively strong 
5h ' sources of gravitational radiation and are expected to be prime sources for the terrestrial network of kilometer-sized 
' interferometric gravitational wave detectors, which will soon be fully operational, or for space-based detectors such 
as the proposed LISA mission In order to detect gravitational radiation and subsequently study the physics of 
these sources it is absolutely necessary to have a prior theoretical knowledge of their dynamics. This is especially 
true because of the method of matched filtering (see |^ for a recent review) that is likely to be employed in order 
to identify true gravitational wave signals "buried" inside the detector's noisy output. The success of this method 
depends on the use of an accurate template of the incoming waveform. 

This paper will focus on the case of extreme mass ratio systems, modelling a massive central object which is a 
spinning (Kerr) black hole while the orbiting body is "light" and compact enough to be considered as a test-particle 
moving in the gravitational field of its companion. There are two important reasons for studying such a model. 

Firstly, due to the extreme mass ratio, the motion of the small mass can be accurately approximated by a geodesic 
trajectory (which is well known [||) and the system's gravitational radiation is well described by first-order black hole 
perturbation theory techniques. The celebrated Teukolsky formalism |^ has proven particularly successful for this 
task. One thus has the opportunity to make a detailed study of a fully relativistic celestial system. For this reason, 
black hole perturbative studies can be used as a test for numerical relativity simulations of two-body systems (and 
vice versa!) 

Secondly, in recent years there has been an accumulation of evidence of the existence of supermassive black holes (of 
mass range 10^ — lO^Mo) in galactic nuclei (including our own Milky Way) jo). It is expected that scattered stellar- 
mass ^ 1 — IOA/q compact objects from the surrounding stellar population will be captured by the central black 
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hole as a result of two-body encounters and interactions with the inhomogeneities of the background gravitational 
potential. The same scenario can of course work equally well for normal stars; however they will soon be tidally 
disrupted as they approach the black hole 0, [|| ||]. 

Once in a bound orbit, the compact object will slowly inspiral towards the central black hole due to the emission 
of gravitational radiation. As the frequency of the emitted waves scales as \/M (where M denotes the central black 
hole's mass), they will potentially lie in the low-frequency band (10^^ — \{)^^Hz) where LISA will have its peak 
sensitivity. The (still uncertain) estimate for the number of such events is around 1/year, or better, out to a distance 
of IGpc and they should be detectable by LISA, with typical signal to noise ratios of 10 — 100, 0, assuming the use 
of some optimal filtering technique, such as matched filtering j^. 

A huge payoff from direct observations of such events is to be expected, provided we have an accurate a priori 
description of the emitted waveform. In principle, for instance, the black hole parameters (masses, spins) can be 
measured to a high accuracy. Similarly, information on the mass-function of compact stellar populations in galactic 
nuclei could be provided. Because the total luminosity of the source depends only on its mass it may be possible 
to work out the distance to the source, which would be very useful to cosmologists. Moreover, one might be able 
to identify the massive object as a Kerr black hole, as opposed to some other, more speculative object (for example 
a boson star [0). This was demonstrated by Ryan [ll| ] who showed in detail how the massive body's multipole 
moments are encrypted in the waveform emitted by an orbiting particle. In the near term, precise numerical results 
in the low-mass-ratio limit will be useful for testing the accuracy of post-Newtonian (PN) derived templates aimed at 
ground based detectors hke LIGO ||l^, [p^ . 

LISA will monitor the last year of inspiral of a compact body into a massive black hole by tracking the phase of the 
emitted waveform. It has been suggested that for astrophysically likely scenarios, drag forces operating on the orbiting 
body due to gas accretingonto the black hole, will operate on a timescale much longer than the radiation reaction 
timescale of the particle [ p^ , jist . Based on requirements that the initial highly eccentric orbit in which the particle 
finds itself as the result of some scattering event should have a small enough periastron so that the radiation reaction 
timescale is shorter than the timescale for a second scattering event at apastron, we expect that the initial periastron 
should be rather close, so that < 20M while the apastron will extend to a distance 10^ — lO^M jl^, 
Newtonian order estimates suggest that although radiation reaction will considerably reduce this enormous initial 
eccentricity during the course of the inspiral, the eccentricity will remain finite and non-negligible when the particle 
enters the strong- field region of interest to this paper (see section VB below). Exactly how much eccentricity remains 
will depend critically on the initial periastron distance and is largely insensitive to the initial apastron distance (and 
thus to the initial eccentricity). 

We can thus argue that for a sufficiently bound orbit the system of the massive black hole and the orbiting compact 
object will evolve under its own spacetime dynamics. This tends to justify our "black hole plus particle" model. 
However, even in this simplified picture there are problems. The particle, in general, will move along a non-equatorial 
eccentric orbit (as the galactic central stellar population is almost spherically symmetric, capture orbits of arbitrary 
inclination are to be expected). The Teukolsky formalism cannot, at present, deal with such orbits, for reasons 
discussed in , in particular the problem of determining the rate of change of the "Carter constant" of the motion 
due to the emission of gravitational waves (much effort, towards this goal, is being focused on building a framework for 
calculating the gravitational self- force acting on the orbiting particle ||l9|). For this reason we restrict our attention to 
equatorial orbits around the central body. In such a case the rate of change of the orbital parameters can be deduced 
by reading the gravitational wave fluxes for the energy and angular momentum at infinity and the black hole horizon. 

There is, however, a factor that cannot be accounted for by the previous flux-balance argument which lies at the 
heart of our approach. As has been observed recently |£l| the gravitational self-force contains a conservative 
piece which is not associated with any radiation emission. Although the effect of this conservative force is negligible 
(scaling as ~ /i^) over short timescales (say, one orbital period), it is conceivable that the same will not be true for 
the accumulated effect after lO"' — 10^ orbits (this is, roughly, the number of orbits that LISA will record). 

Another issue that has risen recently concerns the possible difficulty in defining the notion of adiabaticity and 
averaged flux for generic, i.e. eccentric and non-equatorial, orbits in Kerr spacetime. This is related to the belief 
that generic orbits have no well-defined orbital periods as they show an apparently non-periodic behaviour. For 
that reason, it has been suggested |2^ ] that an "ergodicity" criterion would be more appropriate. However, recent 
work suggests |^ that it is possible, after all, to rigorously define (by means of Hamilton-Jacobi theory) a triplet of 
fundamental frequencies for generic Kerr orbits. Consequently, one may still be able to define adiabaticity for these 
orbits too. 

Serious complications can also arise at the level of free motion, where radiation reaction is neglected. In general, 
the small body will have its own intrinsic spin. In such a case, due to the coupling of the particle's spin with the 
background gravitational field, the motion is no longer geodesic. Although the (specific) spin magnitude is small, i.e. 
S ~ 0{fj,), spin-induced effects could become important over timescales much longer than, say, one orbital period. A 
particularly dramatic possibility is that when the test-particle is allowed to have spin, "chaotic" features may appear 
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in the orbital motion |24|. Presently it is unclear whether chaotic behaviour will be important for extreme mass ratio 
systems likely to be observed by LISA. 

When radiation reaction is "switched-on" in the spinning particle case, one finds, not surprisingly, that the radiative 
fluxes at infinity and the horizon are inadequate for determining the evolution of the orbit. This is, in part, due to 
the fact that there is no known analog of the Carter constant (so there is one less constant of motion available), 
and also due to the existence of additional spin-degrees of freedom. A Newtonian order, weak-field estimation for 
the radiative change of the spin has been worked out by Apostolatos et. al. [P5[ . Some speculations of what could 
happen to circular orbits under strong field conditions can be found in [2^ , p7| . For generic orbits, most likely only 
a self-force calculation will be able to describe the full orbital evolution. 

As we are still far away from dealing with all of these challenges we make two major simplifications for this paper, 
that the orbiting particle has no spin and that it always remains in an equatorial orbit. Although Ryan has shown 
PH , for orbits in the weak field region, that non-equatorial orbits are forced by radiation reaction towards becoming 
retrograde equatorial orbits, the effect is small. The effect remains small even in the strong field region, as was recently 
shown by Hughes [ |T8| . Precisely equatorial orbits, pro- or retrograde, will remain equatorial under radiation reaction. 
Therefore it is reasonable to expect that detectors such as LISA will actually observe signals from particles in near 
equatorial orbits. 

Previous studies have shown that slightly eccentric orbits of particles around Schwarzschild ||2^ and Kerr black 
holes and arbitrarily eccentric orbits around Schwarzschild black holes decrease in their eccentricity until 
shortly before the innermost stable circular orbit (ISCO) when a point is reached after which the eccentricity begins 
to increase. The present work comes as an additional piece to this series of papers. Specifically, we consider equatorial 
eccentric orbits of particles around a Kerr black hole and study their evolution under gravitational radiation reaction. 
This class of orbits is not exactly what we would expect in reality, but it is an important step towards a more realistic 
view of gravitational waves from this type of low-frequency source, because it includes two very important features 
which we know will be present in all or most sources, black hole spin and orbital eccentricity. Our study is, in this 
respect, a useful companion piece to Hughes' discussion of non-equatorial (but circular) orbits 1 18|. Eccentric equatorial 
orbits were first investigated by Shibata who calculated fiuxes and waveforms, without, however, discussing the 
impact of radiation reaction on the orbital motion. Our approach is similar to previous papers investigating eccentric 
orbits around non-spinning black holes , and nearly circular orbits around spinning black holes [^o| and the results 
are qualitatively similar to those of both papers. In addition, we compute gravitational waveforms produced by 
moderate/high eccentricity, strong-field orbits (not discussed in Shibata's study [Q) which we call "zoom- whirl" 
orbits. We find that these waveforms are a very characteristic, though complex, signal that might be important from 
an observational point of view for the planned LISA space antenna. 

In this paper we focus on the final part of the inspiral, when the particle is at small radii, relatively close to the 
last stable bound orbit. In consequence we deal with orbits with moderate eccentricities, between 0.1 and 0.7. In a 
future paper |Q we intend to study the full inspiral, thus expanding our scope to cover orbits with large radii and 
larger eccentricities, on the order of 1. In that paper we plan to present wavetrains and spectra associated with a 
long stretch of the inspiral, covering many orbital periods, along the lines of 

Our results in this paper can be summarised as follows. Moderate eccentricities will be a feature of the signals 
from many inspiralling compact binaries right up to the final plunge. Immediately before plunge there will be an 
eccentricity increasing phase in all cases, particularly noticeable for retrograde orbits. The total amount of eccentricity 
gained in this phase will generally be small, on the order of 10% or less for low-eccentricity (e < 0.1) prograde orbits, 
but perhaps as much as 50% for low-eccentricity retrograde orbits. Orbits with moderate eccentricities will gain much 
less in eccentricity. Where e > 0.3 and the orbit is prograde, zoom-whirl features will be prominent in the waveform 
in the very last stages of the inspiral. Where these orbits are observed from a position away from the polar axis of 
the source there will relatively strong high-frequency component to the signal due to beaming of higher multipoles in 
the radiation in the direction of the orbiting particle's motion. One expects that these signals will present particular 
problems for signal analysis, a situation which may be ameliorated when a positive detection of the source has been 
made during the earlier part of the inspiral when the waveform, though highly eccentric, will be less complex. 



B. Organisation of the paper 



The remainder of this paper is organised as follows. In Part II we discuss the geodesic motion of eccentric equatorial 
orbits (Sections IIA & HB), paying particular attention to the so-called "zoom- whirl" orbits (Section IIC). In those 
sections we define useful orbital parameters such as the semi-latus rectum p and the eccentricity e. Some analytic ap- 
proximations on the orbital periods and number of revolutions for a particle in an orbit close to becoming dynamically 
unstable, are presented in Section IID. Sections IIIA & IIIB of Part III contain a review of the Teukolsky-Sasaki- 
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Nakamura formalism for the calculation of gravitational waveforms and fluxes. In Section IIIC, we give a preliminary 
discussion on the orbital evolution under radiation reaction (definition of adiabaticity, general formulae for the rate 
of change of orbital parameters). Part IV is entirely devoted to analytic results. Section IVA discusses the weak- 
field limit for the orbital parameter's rates of change. In Section IVB we derive an approximate formula relating the 
energy flux to the angular momentum flux, emitted by orbits close to becoming unstable. We subsequently use this 
formula to find strong- field approximate expressions for the rate of change of p and e. In Section IVC we study the 
particularly interesting family of (equatorial) horizon-skimming orbits that can exist around a rapidly rotating black 
hole. The main (numerical) results of this paper are contained in Part V. In Section VA, we sketch the methods 
used in our numerical code and, moreover, give estimates for the various introduced errors. In Section VB we give 
results on the averaged rate of change of the parameters p, e (which determine the evolution of any given orbit). This 
allows us to draw conclusions for the "global" behaviour of bound equatorial orbits under the influence of radiation 
reaction. Section VC contains calculations of waveforms generated from some zoom-whirl orbits. Part V ends with a 
presentation of the spectral content of the radiation emitted at infinity and at the black hole horizon (Section VD). 
Our conclusions are summarised in Part VI, where we also discuss prospects for future work. Tables with samples 
of our numerical data can be found at the end of the paper's main body. Three Appendices are devoted to some 
technical details. Throughout this paper we have adopted geometrised units (c = G = 1). 



II. GEODESIC MOTION 



A. Equations of motion 



We start by considering a test body moving in a Kerr gravitational field. For the moment, we neglect any radiation 
reaction effects and focus on purely geodesic motion. Working in the usual Boyer-Lindquist coordinate frame, the 
equations of motion, specialised for an equatorial orbit, are given by [Q, 

,2j^±(V;)V2, (1) 
r^^^V,^-iaE-L) + f, (2) 



„2 



dt _ , ^ (?- + a )r 



' dr 



V, = -a{aE -L)+ , (3) 



0{r) - tt/2 , (4) 

where T = E{r^ + a^) - La, V,- = - A{r^ + (L - aE)'^), A = - 2Mr + . The two constants of motion E, 
L denote the orbit's specific energy and z-component of angular momentum (for notational simplicity we drop the 
subscript z for the angular momentum). We have prograde (retrograde) orbits according to whether L > (< 0) 
(note that at certain points, where there is no danger of confusion, we shall label retrograde orbits by a negative 
value for the spin parameter a). Moreover, since we shall be discussing bound orbits, < i? < 1. A general bound 
equatorial orbit can be equivalently described Q either by the constants E and L or by a semi-latus rectum p and 
an eccentricity e (with < e < 1). The restriction on the values of p is discussed below. We define these parameters 
in terms of the two turning points of the orbit (vp is the periastron and the apastron, see Fig. ^ for a typical 
illustration) , 

= TT— ' = -r— ■ (5) 



A turning point by definition satisfies Vr{ro) — 0, or explicitly, 

{E"^ - l)r^ + 2Afr2 - [x^ + + 2aEx)r + 2Mx^ = , (6) 

where we have further defined x = L — aE. Writing this polynomial in the form (_E^ — l)(r — rp)(r — ra){r — r^) we 
can immediately write an expression for the energy, 



E 



^2 



,2 



1/2 

(7) 
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Similarly, the third root of m) is found to be 



_ 2A/(1 - e'^)x'^ 
" p2(l_£;2) • 

2F{p,e) 



It then follows that, 



(8) 



(9) 



The explicit forms of the functions N, F and Aj, are given in Appendix A. In this expression, the upper (lower) sign 
corresponds to a prograde (retrograde) orbit. The same convention will be followed throughout the paper. 




FIG. 1. The radial potential Vr (in units of M^**) as a function of r (in units of M) for p — 2.2M, e = 0.5. The black hole 
spin is a = 0.99M. Motion is permissible at the regimes where Vr > 0. It is easy to distinguish the apastron at Va ~ 4.4M and 
the periastron at = 1.47M. The event horizon is at r+ = 1.141M. 



The radial coordinate can be parameterised as 



P 



1 + e cos X 



(10) 



where x is a monotonically varying parameter, running from x = (at r = rp) to x = tt (at r = r^) and finally up to 
X — (back to r = rp ). The radial motion can be separated into two distinct branches, namely, the motion from 
Tp to ra, and the "inverse" motion from ra back to rp again. Integration of (|^) gives 



<(r) 



first branch 



£(r) 

Tj. — t{r) second branch 



where 



iir) 



1 fdr 
V dr 



A 



-{Er^ - ax) 



dr 



We have also denoted as Tr the period of the radial motion. For the 0-motion we similarly write, 

0(r) : 

where 



(/)(r) first branch 

A(p — ip{r) second branch , 



(r) 



1 /dr\"^ 
rf7 j 



X + — {Er'^ ~ ax) 
A^ ' 



dr . 



(11) 



(12) 



(13) 



(14) 



and A(j) is the change of 4> during an interval T^. Both integrands in (|l If) , ( |13| ) are (unphysically) divergent at the 
turning points, an undesirable feature in a numerical calculation. This difficulty can be avoided by choosing x as the 
integration parameter. Using 
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dr esinv, , , 2Mx'^ 



dr p p 

we get 



x^+a^ + 2xaE (3 + e cos x)]^/^ , (15) 



nx) = / dx — \~i/2, : ' (16) 

Jo J{x',p,e)V, (x',p,e) 



Jo J(x',P,e)K- (x ,P,e) 



where 



2Mx^ 

Vr{x,P,e) = + + 2aa;£' (3 + ecosx) , (18) 

P 

2A'Ix 

%ix,Pie) = X + aE —(1 + ecosx), (19) 

9„ '2aMx, , Ep^ 

Vt{x,P,e)^a^E i + e cos x + 7— -^^ , 20 

p (1 + ecosx) 

Jix,P,e) = 1 - —(1 + ecosx) + ^(1 + ecosx)' ■ (21) 
p p'^ 

The integrand quantities in ([l6|), ( |l7| ) are well behaved and, moreover, these equations are valid for both branches of 
the radial motion. The radial period is simply given by = t(27r) = 2i(7r), and similarly, = (t>{2n) — 2(/)(7r). 

A general bound equatorial orbit is the combination of two separable motions: the radial motion which is, strictly 
speaking, periodic (in the sense that the radial coordinate returns to its original value after a certain time interval Tr 
has elapsed) and the azimuthal motion which is not purely periodic (in the sense that the 0-coordinate monotonically 
increases but, nevertheless, the orbit returns to the same configuration after <j) has increased by some value A^). The 
former motion is known in classical mechanics as "libration" while the latter motion is called "rotation" . For such 
a combination of motions, it is generally known that there is a fundamental period (the period of libration) which 
fully describes the motion (see Appendix B for further details). We shall, therefore, call Tr the orbital period. The 
fact that the orbit is periodic in a strict sense will enable us to rigorously define adiabaticity when radiative effects 
are to be included. 

In line with the foregoing discussion, we define the orbital frequency to be 17^ — 2'k jT^. We can similarly refer to 
the frequency of the 0- motion as fi^ = ^(^jTr. The gravitational waves emitted by our systems will have frequencies 
which depend on these orbital frequencies. Below we shall see how they form a spectrum of discrete frequencies 
parameterised by the following wave numbers: which identifies the multipole of the emitted waves {i = 2 for 
quadrupole, for instance), m which runs from — £ to +£ and k which counts the harmonics created by the linear 
composition of the two orbital frequencies. The frequency of the waves emitted by a given harmonic /c of a given 
multipolar contribution m is 

, „ mlS.(h „ 

OJ ^kVtr + ——nr. (22) 

27r 

In the calculational scheme to be outlined below we will evaluate the fluxes of energy and angular momentum which 
are carried by waves of a given frequency (that is, a given multipole and harmonic of the frequency spectrum) and 
sum the fluxes for all frequencies of the discreet spectrum to get the total radiated fluxes of these quantities. 



B. Separatrix curve 

In general eqn. (^ has three distinct real roots. The case with — r^ corresponds to a marginally stable orbit: 
once at the periastron, the particle will enter into a circular orbit of radius rjsbo — — r^ (ISBO stands for Innermost 
Stable Bound Orbit). At this stage the orbit has become unstable, so that a slight inwards "push" will drive the 
particle to catastrophically plunge into the black hole. Therefore, stable bound orbits should satisfy r^ < rp. This 
translates to the inequality, 
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^(l + e)(3-e) <p2 



(23) 



We can imagine a division of the (p, e) plane into regions of stable and unstable orbits. The boundary curve Psie) 
satisfying the equality in (|2^), defines the separatrix of bound orbits. In Fig. |2|we illustrate separatrices for a variety 
of black hole spins. A sample of numerical data used to generate this figure can be found in Table |. As one might 
have anticipated, spinning up the black hole will cause the separatrix curve for prograde (retrograde) orbits to move 
to the left (right) with respect to the Schwarzschild curve Ps{a = 0) = (6 + 2e)M |31|. This behaviour can be seen 
most easily by a slow rotation approximation to (ESf). At leading order we find, 



= (6 + 2e)M =F 8a 



1 



6 + 2e 



1/2 



+ 0(a2) 



(24) 



On the other hand, as can be verified by direct substitution in (23), for extreme rotation (a — M) the prograde 
separatrix becomes Ps{e) — M{1 + e), i.e. for all eccentricities, the periastron "descends" into the black hole "throat" 
at r = M, but is still separated by a finite proper distance from the horizon itself BS]. 




FIG. 2. Separatrices on the (p, e) plane for a variety of black hole spins. From left to right: a/M 
0.5, 0.1, O(dashed), —0.5. As a ^ M the prograde separatrix goes to the limiting value ps M{1 + e). 



0.999,0.99, 



25 



20 



15 



10 



-10 




-25 



-20 



-15 



-10 



10 



15 



20 



25 



FIG. 3. A zoom-whirl orbit with p = 2.35Af, e = 0.9 around an a = 0.99M Kerr black hole. In this figure, the particle has 
performed more than twenty revolutions in less than three orbital periods. The periastron is at rp — 1.237M, located close to 
the hole's event horizon at r-^ = 1.141i\f (denoted by the dashed line). The ISBO radius is risbo = 1.216M. 
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C. Zoom- Whirl orbits 



From the short discussion in the previous Section one can imagine that as the orbit graduaUy approaches the 
separatrix, the particle will spend a considerable amount of its orbital "life" close to the periastron (see Fig. ||) . An 
approximation for as p Ps, derived in the following Section, gives 

T,^-Hp~p,), (25) 

which shows that the period will grow (and eventually diverge) as the separatrix is approached. In that region, 
the particle will trace a quasi-circular path before being reflected back to the apastron. Such behaviour will be 
particularly prominent for high eccentricity orbits: the particle will "zoom in" from its apastron position, and perform 
a certain number of quasi-circular revolutions ( "whirls" ) reaching the periastron (which should have a value close to 
nsbo(e) = ps(e)/(l -|- e)). Finally, the particle will be reflected and "zoom out" towards the apastron again. We shall 
heuristically (but quite descriptively) name these orbits "Zoom- Whirl" orbits. They resemble a set of orbits known in 
the literature as homoclinic orbits [p7| . Zoom- whirl orbits can exist in both Kerr and Schwarzschild geometries, and 
their potential significance for the detection of gravitational waves by space-based instruments was first pointed out 
some years ago by Curt Cutler and Eric PoissonH, who concluded that the small number of whirls in the Schwarzschild 
case made the phenomenon less interesting for spinless central bodies. But as we shall shortly see, they are more 
pronounced in the case of near-extreme Kerr black holes, for prograde orbits. A typical example of such an orbit is 
illustrated in Fig. ||, for the case of a rapidly spinning (a = 0.99M) black hole. 

It is straightforward to calculate the total number of azimuthal revolutions Nr = A(j)/2'K during one orbital period, 



by numerically integrating (16). Results obtained by such a calculation are presented in Fig. |J. In this figure we 
have considered orbits of a given eccentricity (e = 0.9 and e — 0.3) and for a variety of black hole spins. For all 
depicted cases, the smallest value of p resides at the same distance dp from the corresponding separatrix value Ps{e). 
As can be seen, the number of revolutions increases as the separatrix is approached, in agreement with our intuitive 
expectations. In fact, an approximate formula (valid for p pg) derived in Section IID shows that, 

N, ^ - \n{p - p,) . (26) 

We can furthermore deduce that the "whirling" of the particle near the separatrix becomes more pronounced as the 
black hole spin increases. Although for small and moderate spins Nr stays close to the corresponding Schwarzschild 
value, it grows rapidly as a ^ M, basically due to the intense "frame-dragging" induced by the black hole's rotation 
in the very strong field region close to the horizon which can be reached by particles in prograde orbits. The overall 
behaviour can be understood as an extreme example of perihelion advance (as in the celebrated case of the planet 
Mercury) . 



In principle, as (26) suggests, the number of revolutions can be made arbitrarily large irrespective of the black hole 
spin, provided the particle approaches the separatrix sufficiently closely. However, as we discuss in Section IIIC the 
adiabatic assumption upon which our formalism relies breaks down in this regime. Sufficiently close to the separatrix, 
radiation reaction makes a significant correction to the particle's motion in each orbital period. Before long this 
causes the particle to cross the separatrix and plunge into the black hole. These transition/plunging regimes have 
been studied recently by Ori and Thorne for the case of circular equatorial orbits in the Kerr geometry. More 
relevant to the present discussion is the work of O'Shaughnessy and Thorne l39[| which concerns the transition regime 
of zoom- whirl orbits. They show that for the case of an extreme Kerr black hole and eccentricity close to unity, the 
particle may experience more than 20 whirls per orbit before plunging. These have to be added to the number of 
whirls performed during the adiabatic phase of the orbit. 

In a realistic scenario, we should not expect to find (apart from chance cases where the particle enters a near- 
separatrix orbit as a result of its initial scattering) very high eccentricity zoom-whirl orbits, as it is well known that 
the orbit has a general tendency to circularise [Q. However, despite the decrease in eccentricity over the greater 
part of the inspiral, a substantial amount of eccentricity will survive, in many cases, up to the point where the orbit 
is about to plunge. These orbits will probably become zoom-whirl orbits, especially when a rapidly spinning black 
hole is involved and especially for prograde orbits. Keep in mind that many scattered particles will be in highly 
non-equatorial orbits. Zoom-whirl behaviour should also be seen in these cases as there is still a separatrix present, 
close to which the particle can spend a considerable amount of time. 



^The name "Zoom- Whirl" originated with the work of these two at Caltech. It may have been suggested by Kip Thorne. 
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A compact body in a zoom-whirl orbit will spend a considerable fraction of the orbital period in strong field regions 
(it can even travel close to the event horizon if the central black hole is spinning rapidly enough) and hence will 
radiate strongly. Our numerical results together with analytic approximations, reveal that a good fraction of the 
averaged flux is radiated during the motion near the periastron. As the orbit approaches the separatrix it tends to 
radiate as if it was a circular orbit of angular frequency fi^ (see also JsTI for a similar statement in the Schwarzschild 
case). This is clear evidence that most of the radiation is coming from the whirl part of the orbit, during which the 
radius hardly changes and there is a single dominant frequency characterised by the azimuthal (0-dependent) orbital 
period. However, the most important feature of a zoom-whirl orbit is the characteristic form of the gravitational wave 
it emits, which is a series of rapid "quasi-circular" oscillations separated by relatively "quiet" intervals. In Section 
VC below we calculate some waveforms of this type. 

25 I , , , , , , , , 1 

20 - 




' • • ' • • ' • • 1 

23456789 10 
p/M 

30 I , , , , , , , 1 



25 

I 2° - 




I • ' • • ' • • 1 

123456789 
p/M 

FIG. 4. Number of revolutions as a function of the semi-latus rectum p for fixed eccentricity e = 0.9 (top frame) and 
e = 0.3 (bottom frame). The black hole spin is, from right to left, a/M — 0,0.1,0.5,0.99,0.999. Each curve terminates at a 
point located 5p = O.OIA/ away from the respective separatrix value. Evidently, zoom-whirl orbits are expected to be more 
pronounced for rapidly rotating black holes. 



D. Approximations near the separatrix (I) 

Orbits that reside near the separatrix of the (p, e) plane are amenable to analytic approximation, basically due to 
the fact that the turning point rp is close to a local minimum of the radial potential Vr- In this Section we derive 



approximate expressions for and A0. We already know that (see eqn. (17)) 



T. = 2I ^ ^^^^;f> ^ , (27) 



A^ = 2rdx ^^^^ . (28) 

^0 J{x,p,e)Vr {x,P,e) 

We take the "distance" e = p — ps from the separatrix to be small, i.e. e/M ^ 1. Also, we shall exclude small 
eccentricity orbits (more precisely, orbits with e ^ ^/M), or marginally bound orbits (e 1). The former case of 



10 



nearly circular orbits has already been discussed in [Q. In what follows, quantities with an "s" subscript are to be 
evaluated exactly at the separatrix. From (|l8| ) we get, 



where 



K-(X,P, e) - — {1 + 0(e)} [eS + 2exUl - cosx) + e(l - cosx)] 

Ps 



5 = 2p,-(l + e)(3-e) ( — 



(29) 



(30) 



P=Ps 



We see that K- ^ as e ^ and x = (i.e. the periastron "touches" the separatrix). At the same limit, Vt and J 
remain nonzero. We can write then, at leading order in e. 



At{l - cosx) 



[£5 + 26X^(1- cos x)] 



1/2 ' 



(31) 



1/2 



A^{1 - cosx) 



[eS" + 2ea;2(l - cosx)] 



1/2 



We have defined the functions 



My) 



[a'^E,{l + e - eyf - 2aMa:,(l + e - eyf/ps + E^pf\ 
(1 + e - ey)2 [1 - 2M(1 + e - ey)/ps + a2(l + e - eyY/pl] 



(32) 



(33) 



^0(2/) 



[xs + aJ?s - 2Ma:s(l + e - ey)/ps] 
[1 - 2M(1 + e - ey)/ps + a2(l + e - eyf/p^^ ' 



(34) 



with argument y = 1 — cosx- In order to isolate the divergent pieces in the integrals (|3l|), ( ^2|) we split the functions 



A,0(2/) = A,0(o) + Bt,0(y) . 



(35) 



These expressions are just Taylor expansions around the regular point y = Q (with ^ containing the first and all 
higher derivatives of At.^). Not surprisingly, both functions Bt,^{y) take the form 



(36) 



Although we don't write the functions i?t,0(2/) explicitly here (as they don't take a simple form and they are not 
needed in what follows) we have verified that i?t,0(O) 7^ 0. It follows that the contribution to the integrals from 
^t,0(l — cosx) is finite when e,x ^ 0. On the other hand, the contribution from At^^[Q) is found to be divergent at 
the same limit, 



1 



_ i(ea;2)-i/2ln 
[e^ + 2ea;2(l-cosx)]'^' 2 ^ 

Hence at leading order in e (therefore close to the separatrix), 

-(l + e)(3-e)1^/^ 



Tr « A(0) 



eMp, 



In 



eS 



£S'(l + e)(3-e) 



(37) 



(38) 



A0«^0(O) 



(l + e)(3-e) 
eMps 



1/2 



In 



Mepi 



eS'(l + e)(3-e) 



(39) 



The divergence of and At/) at the separatrix is the result of the particle being trapped in an unstable circular 
orbit at the location of the minimum of the radial potential 14 . 
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III. RADIATION REACTION: FORMULATION OF THE PROBLEM 

A. The Teukolsky formalism 

In this paper, we shall employ Teukolsky's formalism Q for the calculation of gravitational fluxes and waveforms. His 
eponymous equation describes the evolution of linearised radiative perturbative fields in a Kerr geometry background. 
In particular, instead of dealing directly with metric perturbations, the Teukolsky formalism considers perturbations 
on the Weyl curvature scalar -04. This quantity is a result of the projection of the Weyl tensor on the null vectors n", 
fh^ which are members of the Newman-Penrose null tetrad [Q, that is ipi = —Caf3'Ysn°'fh^n''fh^. The feature that 
makes this formalism attractive to our problem is that the radiative fluxes (at infinity and at the horizon) as well as 
the two wave polarisations can all be extracted from ■04 • The "master" perturbation equation is separable in 

the Fourier domain by means of a decomposition 

Mt,r,9,cf>) =P~^J2j dcje~'^*+'"''^ ^2SZ{0)R,^^{r) , (40) 

where p — r — ia cos 9. The radial function Rimuj{r) satisfies the Teukolsky equation 

2 d I 1 dR(i„i^ 

The potential V{r) is given by 



- ( - ) - V{r)Ri„.^ = r,,„^ . (41) 



V{r) = — ^- + Siujr + A, (42) 

where K ^ {r'^ + a^)uj — ma and A = Eim + o^lo^ — 2amuj. The angular functions -2Si^{9) are s — —2 spin-weighted 
spheroidal harmonics ||4^ which satisfy the following eigenvalue equation, 

Til 

a^Lu'^ cos^ 9 5 h 4acj cos 9 

sin" 9 



1 


d 






M9 




sin 6 





4m cos 



sm' 



2— -4cot26'-2 + £;^„ 



~2SZ = 0. (43) 



We have adopted the following normalisation for the spheroidal harmonics (hereafter we drop the subscript —2 for 
notational simplicity) , 

r \SZ\^sm9d9 = l. (44) 
Jo 

The source term Timtj present in (|4^) is constructed directly from the particle's energy-momentum tensor and this is 
the point where the particle's motion enters explicitly in the perturbation equation. Its explicit form is given below. 
Let us now return to the radial equation (^l|). A particular solution of this equation can be found in terms of two 
independent solutions R^uj ' ^'imuj '^^ ^^'^ homogeneous equation, 

where W the (constant) Wronskian W[A-^^^R^Z<^, A'^^^R'^^J. The solut ions R^£^^, R^muj "^"^^ chosen such as to 
have, respectively, purely ingoing behaviour at the horizon, and purely outgoing behaviour at infinity. Explicitly, 



nin ^ I A^e for r ^ r+ .^gN 

imuj \ ^SQaut^tur* ^^-lgin^-tu>r- for r ^ +00, ^ ' 



"^irau 1 r^e'"'-' forr^+oo ' ^ ' 
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where k = uj — ma/2Mr+, r_|_ = M + (M^ — a^)^/^ is the outer event horizon, and is the usual tortoise coordinate 
defined by dr^/dr — (r^ + a^)/A. From these expressions we have that W = 2iLoB™. The solution (^ ) describes 
ingoing waves at the horizon and outgoing waves at infinity as it should be required on physical grounds. That is, 

nimu,{r r+) -> 2iujB"' X A2(r') ^ ^Cm^^ \^)(^ I't*) 



The source term Tf^muj is given by |h2| 



r,„„^ - 4 / dndtp-^p-\B'^ + ij^*)e-™'^+-*-l^, (50) 



where 



^2 = -^p'pL^i[p-'Lo{p-'p-'T,,r.)] 

^ p«pA2L_i[p-V^+(p"'p"'A-iT^„)], (51) 



2^/2 



B'* = -^p^pA'j+[p-^J+ip''pT^ 

^ „8-a2t r„-4-2A-l7- /„-2--2 



p''pA'j+[p-''p'A-'L-i{p-'p- Trnn)] ■ (52) 



We have defined the operators 



2V2 



Ls = dg -\ ~ aujsmO + s cot 0, (53) 

sm 



J+ =dr + iK/A . (54) 

(55) 

The quantities T„„, Tfnn,Tfnm are the result of the projection of the particle's energy-momentum tensor T'^'^ on 
the tetrad vectors, i.e. T„„ = T^i^n'^n" etc. The energy-momentum tensor for a particle in an arbitrary orbit 
(t,r(t),0(t),</.(t)) is given by 

2j sm au 

where u'^ = dx'^/dr and Y. — r'^ + cos 9. We obtain for the individual projections |T2[ , 

Tun = f^^6{r - r{t))5{e - 9{t))5{cf> - cf>{t)), (57) 
smt^ 

Tmu = M^'5(r - r{t))Sie - 0(i))<5(</> - 0(i)), (58) 
= M^^(^ - Ki))'5(^ - ^(i))'^('/> - 0(0) , (59) 

with 



Cnn = ^ [E{r^ + «') - «i + ' 



Cmn — 



2^/2T:■ 



:[v}y^ [E{r^ +a^)-aL + T.u'] 



i sin 9 ( ai? — 



L 



C, 



mm — V" / 



i sin 9 ( a_E 



sin^e 



9(0) 



The quantity B(^) represents the effective latitudinal potential, i.e, (Su^)^ — 0{6). Substitution in d5Q) yields 



4/i 



dt / rf6le''^*'™'^(*) 

oo Jo 



+ ^4{p'^£™(p'p"'),r}C^nAp-'r''5(^ - r{tme - 9{t)) 
- \p^A'SZJ+{p''J+{PP~'C^^r^^Sir - r{t))6i9 - 9(1)))} 



where 



Li = d,- 



sm9 



+ auj sin 9 + s cot 9 . 



The ^-integration can be performed directly to give 



where 



Tirauj = M / die 



{Anna + ArfinO + Anmo)'5(?' - r{t)) 



+ {(^mnl + Arnml)6ir - r{t))} ^ + {Arnm2S{r - r{t))} 



e=e{t) 



AnnO — 



-2 



27rA2 



Cr.np-^-p-^L+{p-^L+{p^SZ)}. 



■^rh n 1 



0FA 



CynnP 



V a" 



-asin0(i)5,"--(p-p) 



^ pCmmS 



2 



[ij + *a sin 9{t) [p - p)SZ] 



P ^ pCfnmSZij'-^ + P), 



27r 
1 



P'-^mm'^i 



em- 



Note that all functions of 6* are evaluated at 9 = 9{t). The amplitudes defined in 



'^tmuj 



^) can be written as 



^£miu 
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where 



At] 



dm 



dr 



10} 



dr^ 



r=r{t),e=9{t) 



(72) 



'T-OO 



^ r—r{t),9—9{t) 



dr dr^ 

Up to this point, all expressions listed in this Section are valid for an arbitrary orbit. We now specialise our 
discussion to equatorial orbits by setting 0{t) — -k 12. In this case, I"^^ are functions of r{t) only. As discussed in 
detail in Appendix B, the quantities 

a°°-"{t) ^I°°^"{r{t)) e-""['^(*)-f^**l , (74) 

are periodic functions of time (with period equal to Tj.). Consequently, they can be expanded in a Fourier series 

+00 

a (^) 2^ "fc e - , (75) 

with r^i- = 27r/Ti.. Inverting, we obtain the Fourier coefficients 

a^'" = ^ r dt a^'"{t)e^''''-' . (76) 

Jo 

Using the Fourier series ( |75| ) in ([70|),(71) we arrive at 

+00 

k— — oo 

where tUmk = inQd) + kQr and 



- ^-7%r / dt I-^"{r{t)) ™^(*) . (78) 

Due to symmetries of the Teukolsky equation (^) we have the following conjugation relation, 

^r-™.-. = i-^y^rz ' (80) 

where an overbar denotes complex conjugate. 

We proceed by writing (|78|) as an integral over x. 



r,oo,H A^^r 

^0. 



^ Irak 



dx , ,^V^iL J ^^'cJ (KX)) e-'"^*(>^)-*'"'^('^) . (81) 



2iuj„,kB'" Jo j(^) f//^(x) 



As in the case of eqns. (p^, (p^), this expression is well-behaved at the orbital turning points. Moreover, noting that 
the x-dependence of the integrand in (|8l| ) appears in terms of the form cosx (in terms with r(x))) and sinx (in terms 
with u' ) we can write, 



H_ fi^r r, ^t(x) 

Zj, 



+ ^/~fr-i(Kx))e-^-'"'=*(^H*™^(^)l . (82) 
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The subscripts (±) mean "sinx ±sinx" in the functions 1'^^ ■ The numerical calculation of the amplitudes 
Zj^^ is the "backbone" of our radiation reaction code (see Section V) . Finally, we express the A and C amplitudes 
(pq),®, in terms of x,P, e, 



and 



Vt4(x,P, e) 



2V2pWtix,P,e) 
x^J{x.P,e) 



p^E — ax{l + ecosx)^ + epsinxV^^^'^ix^Pi 
(1 + ecosx) P^E — ax{l + ecosx)^ + epsmxVy^{x,P, e) 



2pWtix,P,e) 



(1 + ecosx) 



2a2u'^ + {ia{auj - m) - 41/}^^ + 2u + ; 



de 



'-{^/2) + {au-m)SZ{Ti/2) 



(83) 
(84) 
(85) 



(86) 



Ammoiu) 



1 _ Qauil lty\ 

, \-2ia^{aL0 - m)u^ + a(acj - m){6iM + a(a^ - m)}u^ 



2^u2(i _ 2Afw + a2u2)2 
4ia(aw — m)u'^ + 2uj{iM + a{auj — m)}u^ — 2iuju + uj'^] , 



(87) 



Afnnl (u) 



Ctrl') 



/^u(l - 2Mu + a2u2) 



neat 

de 



-{TTl2) + {auj~m)SZ{^m 



Afiiml{u) = - 



TT u^(l — 2Mu + a^M-^) ■' 



(89) 



^mm2(w) 



1 CmrhSil^i'n /2) 

/2tt u2 



(90) 



Anno{u) 



TT (1 - 2Afw + a2u2)2 



-2ia 



dS: 



V de 



''-{T:/2) + {aoj-m)SZ{^m]u 



pp, oato 

+ (V2) + 1[au m)^(V2) + {(ac. - m)2 - 2}5,«-(V2) 



961 



de 



(91) 



where m(x,p, e) = (1 + ecosx)/p. By means of (|77| ) one can obtain the following expressions for ■04 at infinity and on 
the horizon, 



where 



ij}i{t,r,e,(t)) 



I 7/,oo 
Vimk 



12-K 



Eimk ^Z,k for r -> +00 , 

yoo,H natLJrnk /l^ „ — ii^„, fc ( t — r * ) +im 



(92) 



(93) 



Once the Weyl scalar V'4 is known, we can immediately relate it to the two polarisation components , /ix of the 
transverse-traceless metric perturbation as r ^ cx3 M , 
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1 fd^h+ d^h 



It follows from (Q), (^) that ft,+.x(i, f", ^, 4>) are given by (here the coordinates t, r, 9, (j) are referred to the observation 
point), 

9 yH qa^mk f o\ 

h+-lK^-y ±I^^Irn^^-^^r..(t-r,)+^m^ (95) 

Note that the gravitational waveform is exclusively radiated at harmonics of the two orbital frequencies 17^, ^4,- The 
gravitational wave energy and angular momentum flux at infinity can be found in terms of the Landau-Lifschitz 
pseudotensor Gsf, 



GW 



dE 
~dt 

dL 

/ GW 

We define the averaged (over one orbital period) fluxes to be (C = E,L), 

GW 




r^dn , (96) 

(97) 



With the help of (^, (|97|), (|95|) we arrive at |4| 

^GW = E ' 

- 2^ 4^^3 • (100) 

The calculation of the respective fluxes at the black hole horizon is a more complicated issue as it is not possible to 
use expressions such as (p6|), (p7|). Despite this difficulty, Teukolsky and Press |44j were able to derive formulae for 
the horizon fluxes using the approach of Hawking and Hartle |^ who studied the deformation of the hole's event 
horizon under the influence of infalling radiation. These formulae are, 

^GW = 2^ atmk- — — , (101) 

e,m,k ^^rnk 

tH „ ™l^imfeP ^^(\^\ 
^GW = Otimk—. 3 , (1U2) 



e,m,k '^'^mk 



where 



256(2Mr+)5p^fe(p^, + 4£^)(p^, + 16e^)a;^, 



with e = VAf 2 - a2/4Mr+ and 



Ctoife = [(A + 2)^ + Aamujmk - 4a^w^fe](A^ + 36aTOWmfe - iQa^uj"^^) 

+ (2A + 3)(96a2w^j^ - ASamuj^k) + 144w2„JM2 ~ a^) , (104) 



is the so-called Starobinsky constant. Note that eqns. (|99[)-(102) have to be divided by /i in order to convert them 
to fluxes of specific energy and angular momentum. Moreover, we can exploit the conjugation relations ( pO[ ) in the 
numerical calculation of the amplitudes Z(,mk and reduce by one half the required computational time. 
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B. The Sasaki-Nakamura equation 



From cqn. ( p2| ) it is obvious th at in order to calculate the amplitudes Z'^j^j^ , which will give us the gravitational 
waveform and fluxes (p5|), (p9|)- ( |102| ), we need to evaluate the quantity B™. In principle, one could numerically 
integrate the Teukolsky equation ([11[) from the horizon out to "infinity" and extract the amplitudes B™'°^^. But this 
is a poor strategy, since the effective potential V{r) is long-ranged and the term drops off towards infinity much 
faster than the B°^^ term and can only be extracted with very low accuracy |Q . A way to circumvent this difficulty 
is to integrate, instead, the Sasaki-Nakamura equation Il3], 



(fx , ,dX 







(105) 



The "potentials" F{r),U{r) are given in Appendix C. The solutions of this equation are related to the solutions of 
the Teukolsky equation via the transformation. 



A 



(r2 -f a2)i/2 ~ Ad? 



AX. 



(r2 +a2)i/2 



(106) 



The functions a(r), /?(r) are also given in Appendix C. The key property of ( |105| ) is that it encompasses a short-range 
potential. This can be demonstrated more easily if we shift to the function. 



Y{r) = r;i/2(r)A:(r) . 



Then, eqn. (105) transforms into the Schrodinger-type equation. 



drl 



with the effective potential, 



-(?7,,)2 + 2Mr/, 



+ a 



The functions F(r), U{r) have the following behaviour at infinity and at the horizon. 



Fir) 



+ O{r-r+) 



for r 
for r 



(107) 



(108) 



(109) 



(110) 



TT( \ / -A:2 + 0(r - r+) for r ^ r+ , , 

^ \ -cj2 + [A + 2(1 + amtu - c?uJ^) - luoci/c^] + 0{r-^), for r ^ +oo . ^^^^> 



It follows that. 



^, -s 1 (jj2 — 2 -I- amcj — a2(jj2 _ j^d/co] + ©(r, lnr») for — > +cx) l^^o\ 

U2+0(e-.) forn^-oo, ^''^> 

where c = (r+ — r_)/2M is a positive constant. From this expression is obvious that Q is short-ranged. Consequently, 
equation (105) admits a solution ("in" mode) which is purely ingoing at the horizon and a mixture of ingoing/outgoing 
waves at infinity: 



X'" ^ ] A e for r ^ r+ ('['[X] 

' A'"e-*"'^* -I- A°"'e*"''- for r +oo, ^ ' 



Another useful independent solution to (105) is the "up" mode. 



^up _^ ] Ding-*''-'-* -I- D°"te*'''''* for r r+ 



The relation between the asymptotic amphtudes appearing in (ttq) and (113) can be deduced from (UOq) 
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B'" = -^A'^ (115) 

^out ^_l^^out ^ (116) 
Co 



where the constant cq is given in Appendix C. Hence, we can simply integrate equation (105) instead of ([4l|) and 



easily identify the ingoing and ou tgoin g waves and evaluate their respective amplitudes. We can then simply find the 



desired amplitudes 5'"/°"* from (116). Similarly, knowledge of the wavefunction X{r) and its derivative at a given 



point immediately leads to the Teukolsky radial function R{r) and its derivative via the rule (106). In conclusion 



all the quantities (apart from the spheroidal harmonics) required for the calculation of t he gravitational flux and 



waveform, can be obtained by numerical integration of the Sasaki-Nakamura equation (105) 



C. Orbital evolution: adiabaticity and flux balance 

Due to the emission of gravitational radiation the orbit of a particle around a black hole will slowly evolve in time 
and the orbital constants E, L (or equivalently p, e) will no longer be conserved. Radiation reaction effects become 
noticeable on a timescale that scales as ~ M^//i, i.e. they are always tiny in a timescale ~ 0{M), provided the 
system's mass ratio is sufficiently small. We can define as the radiation reaction timescale, 

Trr = min[rp,Te] , (117) 

where Te = e/|e| and Tp = p/\p\ are the radiative timescales for p and e respectively (approximate expressions for 
these timescales are given in the following Section). We will then say that an orbit evolves adiabatically if 

Trr » T, . (118) 

In other words, it is a good approximation to assume the motion of the particle to be geodesic, as long as we are 
interested in timescales much shorter than Trr. On the other hand, by making this simplification we "freeze" the 
evolution of the orbit, as if there was no radiation reaction. That is, within the adiabatic approximation, we cannot 
know the exact evolution of the functions E{t),L{t) (or of p{t),e{t)). It is still possible, however, to calculate an 
averaged rate of change of such quantities. This can be done by assuming the following "flux-balance" relation 

C = -Cgw = -(C'gw + C'gw) I (119) 

where C — E, L. We have separately denoted the gravitational wave fluxes at infinity and down to the horizon by 
^GWj ^gw respectively. The overdot symbol stands for the averaged (over one orbital period) rate, see Section IIIB. 
We can equally well describe an orbit by means of the parameters (p, e) and calculate the relevant averaged rates of 
change of those quantities. Since E = E{p, e) and L — L{p, e) we have that (commas denote partial derivatives), 

E = E^pp + E^ee , (120) 

L^L^pP + L^ee. (121) 

These relations can be inverted to obtain, 

p = H-\-E^^L + L^eE) , (122) 

e = H-\E,pL - L^pE) , (123) 

with H = E^pL^e — E^eL^p. Eventually, all partial derivatives of E and L can be found in terms of the corresponding 
partial derivatives of the functions F, N and Ax which are given explicitly in Appendix A. However, the resulting 
formulae are quite messy so we do not present them here. 

Note that although the formalism adopted in our analysis offers only a "local" information on the radiative orbital 
evolution, it can be further manipulated in order to obtain additional information. A recipe for "evolving" orbits under 
radiation reaction, using the known averaged rates of change of the relevant orbital constants, was given recently by 
Hughes I^J] in the context of circular non-equatorial orbits. In effect, one is able to construct a series of "snapshots" 
of the radiation-induced inspiral, and make predictions of the evolution of the emitted waveform close to the point 
where the orbit becomes unstable (that is until the adiabatic condition no longer holds). 

As we have already mentioned, adiabaticity will eventually break down near the separatrix, no matter how small 



the mass ratio is. This can be immediately seen from (118) and recalling that ^ oo at the separatrix, as predicted 
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by (38). For an order- of- m agnitude estimation, we can use the quadrupole approximation for the fluxes (see next 
Section) and translate (118) into a constraint on the mass ratio, 



/i 5 
M 1287r \M 



in (' 



P 



5/2 



1± 



a ( M 



M \ p 
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(124) 



The function /3(e) is defined in the following Section. Equation (124) is accurate to leading order in M/p and a/M, 
and to derive it we have used the corresponding order expression for the orbital period 
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The mass-ratio constraint (124) is automatically satisfied as the black hole perturbation scheme we employ requires 
H/M < 1. 

On the other hand, in the strong-field regime near the separatrix we find (using results derived in Section IVB), 
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where (5 is a combination of A((0), j40(O), -B^/e, L and is of order unity. As we discuss in Section IVB, the 
quantity H also becomes zero when e ^ (unless a ^ M, in which case it remains finite). This is clearly the most 
severe restriction for the mass ratio. Fortunately, the real astrophysical systems we are trying to model are typically 
characterised by a mass ratio n/M ^ 10~^. Therefore we can approach the separatrix closely, probably to the point 
where the physical body would begin its plunge into the black hole, in the cases which interest us. 



IV. ANALYTICAL RESULTS 



A. Weak-field approximations 



Orbits with p ^ M are well described by weak-field approximate results. In particular, the energy and angular 
momentum fluxes should be given with sufflcient accuracy by the quadrupole-order formulae as given in p^ , [ p2[ . 
However, these authors make use of a different set of orbital parameters. For example, Ryan's semi-major axis a and 
eccentricity e [E8[ are related to our parameters via the transformation. 
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Note that at a Newtonian level the two sets are consistent to each other. Rewriting Ryan's fluxes in terms of our 
parameters we obtain : 
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The formulae (|129| ) , ( 130 ) can be utilised for order of magnitude estimations even in the strong field regime though 
becoming increasingly inaccurate with decreasing p/M (this has been verified by comparing them to the fully numerical 
results ). We can now estimate the timescales of radiative evolution for p, e. For p ^ M the energy and angular 
momentum, at leading order in AI/p and a/M, are given by, 
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Accordingly, eqns. (122), (|123|) become 
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The equations above demonstrate the well known fact |4^] that, in the weak-field regime, eccentric orbits tend to 
circularise under radiation reaction (while they slowly shrink towards the central body). For the associated timescales 
one finds. 
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According to this equation, in the weak-field regime the eccentricity decays faster than the size of the orbit |^. The 
leading order spin term furthermore implies that this behaviour is more pronounced for retrograde orbits. 



B. Approximations near the separatrix (II) 



The previous Section discussed results which are already familiar from the existing literature ||3l[] . We now present 
new results regarding strong-field orbits which reside near the separatrix. 

The analysis of Section IIIA has shown that the gravitational wave spectrum will essentially contain harmonics 
of r^r, ri^. We can use the approximate expressions (^8|),(|39|) to deduce that for orbits near the separatrix, i.e. 
p ~ Ps{e) — e <ti M (and as long as e is not close to zero and is not unity). 
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Hence, for e ^ we have fi^ — > and in effect, the spectrum becomes almost "circular": 



Furthermore, by substitution in (|99|), (IOC) we get 
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We conclude that orbits near the separatrix radiate energy and angular momentum at rates so that the ratio 
Egw/Lqw is almost equal to the respective ratio of a circular orbit with the same fi^. The effective radius of 
this fiducial orbit is 
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For example, for the prograde orbit p = 2.11M, e = 0.7 we find VeS = 1-88M > Vp 
p = 2.35M, e = 0.9 we find reS = 3.90M > = 1.24M (for both cases we have taken a 



(148) 

1.24M while for the orbit 
0.99M ). 

Equation (147) suggests that particles in zoom- whirl orbits lose most of their energy and angular momentum while 
they revolve near the periastron, which is what we would intuitively expect. 

We next discuss approximations for p and e near the separatrix. Unfortunately, the lack of a simple analytic 
expression for Ps{e) makes such a task difficult, and the resulting formulae are quite cumbersome with little analytic 
transparen cy. Neverthe less , w e ca n follow a much simpler path and still gain some significant insight. For p fa ps and 
using eq. (147), eqns. (|122D, (123) become, 
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By direct substitution of (P3[), it can be shown that the function H{p,e) becomes exactly z ero a t the s eparatrix. On 
the other hand, and as long as a 7^ M, one can verify numerically that the numerators in (14£), ( |150| ) remain finite 
near and at the separatrix. It follows that for non-extreme Kerr holes both p and e diverge at the location of the 
separatrix. This pathological behaviour signals the breakdown of the adiabaticity assumption upon which our method 
stands. A proper discussion of this transition regime should take into account the rapid radiative evolution of the 
orbit (which now varies i n a t imesc ale c omparable to the orbital period). 
Moving on, we divide (150) into (149) to get, 
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Exploring the numerical value of this quantity for numerous very near-separatrix orbits and black hole spins a < M, 
we have found it to be always negative and finite. This means that e and p have opposite signs near the separatrix. 
Since the latter is always negative (the orbit always shrinks) we conclude that ver y c lose to the separatrix e > 0, 
i.e. the orbit gains eccentricity, (as previously found, in less general cases, in |Q, |^). Since weak field 
orbits always lose eccentricity, there must be a critical curve Pcritifi) on the (p, e) plane at which e = 0. As eqn. 
( |l5lD is formally accurate (within the constraints imposed by the adiabaticity condition) not only at the separatrix 
but also in its vicinity, we can actually study the be havio ur of e in a thin zone near the separatrix. For a given small 
or moderate black hole spin, we find that the ratio (151) is again negative. However, for high eccentricities e w 1 we 
initially get a positive value which gradually passes from zero and becomes negative as p — > pc. With increasing a/M 
we observe the same behaviour at even lower eccentricities , pr ovided we are considering prograde orbits. The opposite 
behaviour is observed for retrograde orbits. For a ~ M, (151) becomes negative only very close to the separatrix for 
all eccentricities. These results suggest that, at least for e « 1, the critical curve Pcrit{e) is located close to Ps{e) 
(this has been shown to be true for a = |31|), and that (for prograde orbits) Pcrit{e) — Ps{e) — > as a ^ M (which 
resembles the situation for nearly circular equatorial Kerr orbits pO|). All of our (semi) analytic predictions are fully 
supported by the numerical results that are presented in Section VB. 

Prograde orbits near the separatrix of an extreme Kerr hole are discussed separately, and fully analytically, in the 
following Section. Here we should emphasise once again that all the approximations presented in this Section are 
valid provided e ^ e/M. This excludes nearly circular orbits, which have been explored in detail in | |30[ |. 

We can now write approximate expressions for the timescales Te,Tp for an orbit close to the separatrix. 
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For example, for a a = 0.99Af Kerr hole, wc have T^/Tp = 0.81 for p = 1.7M, e = 0.3 while for p = 2.11M, e = 0.7 we 
get Te/Tp — 4.5 (for both orbits, eqn. (152) is a good approximation). This situation is typical for non-extreme holes. 
As we move along the separatrix keeping a fixed distance from it, the ratio T^/Tp tends to increase (and beco me lar ger 
than unity) with increasing eccentricity. In comparison, the corresponding weak-field timescales ratio, eqn. ( |143| ), is 
relatively insensitive to variations of eccentricity. 



C. Horizon-skimming orbits 

A particularly interesting class of prograde strong-field orbits are those that potentially "graze" the black hole 
horizon. These orbits can only exist provided the black hole is near extremally rotating, a w M (this can be deduced 
from Fig. H). Circular, non-equatorial horizon-skimming orbits were first studied by Wilkins ]48| and more recently 
by Hughes |49| ]. Here, on the other hand, we discuss equatorial horizon-skimming orbits of arbitrary (but not equal 
to unity or close to zero) eccentricity around an extreme Kerr black hole. 

As the separatrix for these orbits takes the very simple form Ps{e) — M{1 + e) we can duplicate the analysis of the 
previous section following a purely analytical path. Expanding (^ around p — Ps — A/(l -t- e) we find that 
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We then get for the energy and angular momentum. 
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L = 2A/^^ + /(e)^ + 0{e^) . (156) 
The explicit form of the functions g{e), /(e) is not required for the following analysis. We use these equations (together 



with (147), noting that = 1/2M on the separatrix for the orbits under discussion) to obtain, 
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E,pL - L,pE ^ [2Mgie) - fie)]E . 

Here, unlike the non-extreme case, the function H remains finite as e — *■ 0. The above formulae, as well as the 
following ones, have a fractional error 0{t/eM). Hence for horizon skimming orbits, 

iA/(l-he)i/2(3-e)3/2£;, (160) 

e^\{l + eY/\?,-ef'^E . (161) 

We see that both rates are finite all the way down to the separatrix, unlike the a < M case. However, the adiabaticity 



condition (118) is still invalidated at p = 



More interesting is the behaviour of the ratio of the rates (16C), (161), 

- = i--HO(e/eM). (162) 
p M 

This is always positive, which means that the e > region that exists for a < M shrinks to zero for extreme Kerr 
black holes. In other words, the critical curve Pcrit{e) has the same value of the Boyer-Lindquist coordinate as the 
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separatrix itself. This conclusion completes the discussion of the previous subsection. For the ratio of the respective 
timescales we get, 

| = ^ + a(e/eM). (163) 

As ( |163D predicts, the eccentricity always evolves more rapidly than the senii-latus rectum, which again is contrary 
to the situation with weak-field orbits. 

V. NUMERICAL RESULTS 

A. Method and error estimates 

Numerical solution of the Tcukolsky equation or related equations has been a minor industry for nearly thirty years, 
since the pioneering work of Press and Teukolsky , Detweiler and Sasaki and Nakamura . Our method is 
based on a numerical algorithm outlined in |]5^ and employs subroutines found in . It involves the use of Bulirsch- 
Stoer integration to solve the Sasaki-Nakamura equation. Our code is a direct descendent of the codes used in ]3l| ] 
and |^o|| , since we deal with both arbitrarily eccentric orbits and black holes of non-zero spin. Romberg integration is 
used both to integrate Eqns. (|l^) and (|l^) and to integrate Eq n. (^ ). To calculate the spheroidal harmonic functions 
^im "spectral decomposition" method described in |18[. The reliability of these methods in general is well 

known. 



We were able to check our numerical results for circular equatorial orbits with those of |18|, where our agreement 
was good to 5 or 6 significant digits, and with the codes on which this code was based Q, |30|, for eccentric orbits in 
Schwarzschild and for nearly circular orbits in Kerr, and again our agreement was good to 4 to 6 si gnif icant digits. A 
similarly good agreement was achieved by comparison with the published results of Tanaka et. al. |pO| for equatorial 
eccentric orbits in Schwarzschild and with those of Shibat a |51|| for circular equatorial orbits in Kerr. We were also 
able to compare our results with those given by Shibata |3^ for equatorial eccentric orbits in Kerr. In this case, 
however, we did find some disagreement of about ^ 1%. The cause of this disagreement is not apparent, while it stays 
roughly at the same level for moderate and high eccentricities. The disagreement does not seem to be due to the 
problems of maintaining accuracy with the long runtimes and large number of harmonics required for moderate/large 
eccentricities. The error introduced by the truncation of the i, k sums in the flux calculation does not seem to be the 
source of the disagreement. We cannot say at present which code might be at fault. 

Finally, we have also been able to compare our code with some results for circular orbits in Kerr from , and again 
our agreement is good to several significant figures. Similarly, comparison with post Newtonian results for eccentric 



Kerr orbits (as found, for example, in [2^ ) reveals good agreement in the weak-field regime. In Table III we compare 
some sample results. 

In view of the lack of any check for strong field orbits with high/moderate eccentricities and high spins, it is of 
obvious importance that we present some estimate of the likely error in our numerical results. The main sources of 
numerical error in our code are as follows: 

1. Inaccuracy in the Bulirsch-Stoer integration routines, from fs^ . We set the relative accuracy parameter eps, 
which governs the convergence of the final result, at 10~^. 

2. Inaccuracies in the Romberg integrator, also from |^^. We set eps — 10^^ for the routines which integrated 
Eqs. ( p^ ) and (p^). This parameter governs the level of convergence which the routine demands in the final result, 
before it stops iterating. However, in the case of the Romberg routine which governed the main program loop, i.e. the 
integration of Eq. (|8^), we typically set eps = 10~^ in many cases in order to achieve large savings in computing times. 
In a few runs designed to produce data for illustration of waveforms only (not numerical data on flux quantities), we 
used eps — 10~^. 



3. Our method requires that we calculate the quantity S'" (see Eqn. {p2\) above). To do this we integrated the 



Sasaki-Nakamura equation ( 105 ) out to r = 100 /ujmk and then successively doubled the limit of integration, until our 
Richardson cxtrapolator told us we had achieved convergence to a relative accuracy of 10~^. 

4. Our method for calculating spheroidal harmonic functions involved writing them as an expansion of an 
infinite set of spherical harmonic functions. Fortunately this expansion can be truncated at 30 terms and remain very 
precise in most cases, but for high black hole spins, a and high angular frequencies, ujmk we were obliged to use 40 
terms to avoid truncation errors causing small high frequency ripples in the wave forms. However, in our numerical 
results of flux rate and orbital evolution this source of error appears to be considerably less than 10~^. 

5. In principle our calculation of fluxes must be summed over an infinite number of harmonics in each of the integers 
£,m and k. In practice truncating these sums for the £ and m harmonics was not difficult. Fluxes for a sequence 
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of these harmonics usually monotonically decrease after a peak at some value of I and m and so we demanded that 
the loop through these variables halt once fluxes went below a factor of 10~^ times the peak contribution. However 
the spectrum of fluxes in the k harmonic was much more complicated, typically involving several peaks (see figures 
in Section VD) before finally monotonically decreasing after a number of peaks which increased for increasingly 
eccentric orbits. We examined spectra in k by hand to confirm the machine's results and experimented widely to 
convince ourselves that we had caught all significant contributions to the total flux from different frequencies. 

Clearly there are several significant independent sources of error, so that we can only offer our best judgement of 
the total relative error in our code in those cases where we have no independent check on our results. As we have 
addressed every systematic source of error that we encormtered, and as we are confident that the code is running 
correctly in all of the cases dealt with in this paper, we estimate the relative error for numerical results quoted in this 
paper as no greater than 10^"^ — 10^"', in the case of fiuxes, E and L, and no greater than 10~^ — 10"'^ for quantities 
such as e and p because cancellations between terms when converting flux numbers into orbital evolution quantities 
tends to increase the size of the relative errors. This is especially true near the critical point where the rate of change 
of the eccentricity becomes zero, due to the complete cancellation of these terms. As an illustration of this, we will 
note in passing that the mysterious "bump" seen in Fig. 2 of Ref. |^ turns out to be due to a rare case where 
flux errors which appear insignificant themselves are greatly magnified when the numerical flux data is combined to 
produce e and p. 

It is useful, in this context, to mention that in comparing the results from our code to the code in the flux 
quantities for radiation emitted toward infinity agree to about 1 part in 10^®, the flux quantities for radiation towards 
the horizon agree to about 1 part in 10^'* and the position of the critical curve, as calculated by the two codes, can 
disagree by about 1%. This suggests that the only way in which our numerical error is large enough to make a visible 
difference in our figures would be as a slight change in the position of the critical curve in Fig. || (the retrograde case). 

B. Backreaction on the orbit 

In this Section we present numerical results on the evolution of bound equatorial orbits in terms of p and e. This 
pair of parameters is preferable to the equivalent set E, L, because of their clearer geometrical meaning. We have 
calculated the averaged rates p,e for a number of prograde and retrograde orbits and for two different black hole spins: 
a = 0.5M, and a = 0.99M. A representative part of our numerical results is presented in Fig. ^ Each individual orbit 
is represented as a point on the (p, e) plane. At each one of those points we have attached a vector with components 
{M/^){p, Me) that indicates the direction at which the orbit adiabatically evolves under radiation reaction. Moreover, 
all orbits shown are chosen so as to be strongly adiabatic for the typical mass ratio 10~^ (the most severe constraint 
is n/M ^ 10^^ ). These figures (together with the analytic approximations of Sections IVA, IVB and IVC) lead to 
the following conclusions : 

• The semi-latus rectum p always decreases (the orbit is shrinking). In fact, \p\ grows monotonically, and finally 
diverges, as the separatrix is reached. This divergence is an artificial feature of our formalism, associated to the 
breakdown of the adiabatic approximation. 

• The eccentricity e shows a more complicated behaviour. For sufficiently large p, we always find e < 0. However, 
as the orbit approaches the separatrix, e changes sign and becomes positive, i.e near the separatrix, eccentricity 
increases. As in the case of p, e will also diverge at the separatrix, due to the failure of adiabaticity. 

• As the black hole spins up, in the case of prograde orbits the critical radius after which e > moves closer to 
the separatrix (in coordinate terms) for a given e. The same is true for a fixed black hole spin, but for increasing 
e. For retrograde orbits, spinning up the black hole tends to move the e = curve away from the separatrix. 

• In a sense, the increasing eccentricity regime is a precursor of orbital instability and plunging. This is hinted by 
the proximity of the critical curve Pcrit, where e flips sign and becomes positive, to the separatrix curve Ps(e) 
which is the boundary between stable and unstable bound orbits. Qualitively speaking, at this stage of the 
inspiral, the radial potential Vr is quite "flat" and as a consequence the particle has more room to move radially, 
even as it continues to "sink" towards the bottom of the potential well (the behaviour which is responsible for 
the characteristic "circularizing" tendency). 

These results agree with, and at the same time generalise previous results concerning bound orbits (of arbitrary e) 
around Schwarzschild black holes |3l) and slightly eccentric equatorial orbits around Kerr black holes |3^] . As we have 
discussed in Section VC, some of the conclusions above must be modified when the black hole is extreme (a — M). 
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In Table VII we give a sample of our numerical data, for the energy and angular momentum fluxes as well as for 
p, e, for some of the orbits presented in Fig. |[ As we have discussed, we believe that these numbers have fractional 
accuracy at least 10^'^. 
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FIG. 5. The evolution of a family of eccentric equatorial orbits, illustrated on the (p, e) plane. The black hole spin is 
a = 0.99A'/ (top graph, prograde orbits; bottom graph, retrograde orbits) and a = 0.5M (graph in the middle; prograde orbits). 
The dashed curves represents the separatrix of stable orbits, while the solod curves represents the critical e = curve. Each 
orbit corresponds to a point in the graph, and its (adiabatic) radiative evolution is represented by a vector with components 
{M/fi){p, Me). Solid and dashed arrows represent the orbital evolution respectively with and without including the fluxes at the 
black hole horizon (the difference between these arrows is visible only in the a — 0.99A/ case). When the black hole is rapidly 
spinning the horizon flux effectively represents gain of energy for the orbiting particle, an effect attributed to superradiance, 
see discussion in the main text for more details. As a result, the inspiral of the body is stalled and the critical curve is 
slightly pushed outwards. This is the reason for the strong misalignment between the solid and dashed vectors at the point 
p = 1.9, e = 0.5 of the a = 0.99A/ plane. Note the much more pronounced orbital evolution for the prograde a — 0.99Af case 
(a consequence of the particle's motion in very strong field regions) and the approach (diverge) of the relative positions of the 
separatrix and the critical curve between the a = 0.5M case and the prograde (retrograde) a — 0.99A'/ case. 



Another important result concerns the significance of the horizon fluxes on the evolution of orbits with relatively 
small periastrii. Specifically, we have encountered very-strong field orbits for which I^^cvfI However, the 

most intriguing property of the horizon fluxes is that they assist, in a sense, the orbiting body. This is most easily 
illustrated by plotting the evolution of the set of orbits of the top graph of Fig. ||, without including the horizon 
flux (represented by dashed arrows in Fig. ^; the solid arrows represent the total rates). For very strong field orbits, 
when the horizon fluxes are taken into account, the shrinking of the orbit is noticeably stalled. Typically, when 
Eq^y is non-negligible it also happens that it represents energy gain instead of energy loss (in other words the fluxes 
EQ^r,EQ'^r have opposite signs). The orbiting particle is effectively draining energy from the black hole itself. This 
is just a manifestation of the so-called superradiance phenomenon, well known in black hole physics [Q, |5^: waves 
scattered within the black hole's ergoregion and having frequencies (as measured at infinity) that lie in the interval 
< uj < 77ia;_|_, effectively appear (for a distant observer) as emerging from the horizon, and amplified at the expense 
of the hole's rotational energy. The outgoing reflected waves "push" the particle outwards, and this interaction is 
manifested as a gain of orbital energy and angular momentum. Our result can be easily understood if we recall that the 
ergoregion is growing for increasing black hole spin. At the same time, because the boundary of instability moves in to 
lower radii with increasing spin, a particle can enter regions with much stronger flelds and therefore emit a substantial 
amount of radiation towards the ergoregion (see Fig. 0). Hence we flnd a signiflcant negative (superradiant) horizon 
flux. An alternative way of viewing this phenomenon is as an exchange of energy and angular momentum via tidal 
coupling analogous to tidal friction in the Earth- Moon system (and elsewhere). For an exposition of this intuitively 
instructive viewpoint see , and references therein. 

It was recognised long ago that if the superradiance effect ever became large enough a floating orbit would 
result when it balanced the energy loss due to radiation emitted towards infinity. Our results confirm earlier work 
in collaboration with Scott Hughes, suggesting that even for orbits very close to the horizon of very rapidly rotating 
black holes, the gain in energy from superradiance is only 10% of the energy lost by the system as a whole. For further 
details, see fsll. 
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The results presented so far in this Section, ahhough very insightful, are still incomplete in the sense that they do 
not describe the radiative evolution of a single given orbit on the {p, e) plane. Instead, they provide local information 
on the evolution of an orbit at a given point. An effort to "paste" together a sequence of such points, in order to 
follow the full orbit is currently underway |]33|] . Meanwhile, we can use certain approximations to foresee what the 
full inspiral trajectory will look like. As a starting point we use the leading quadrupole-order expressions for p, e, i.e. 



setting a = in eqns (|138[ ), (|139| ) and derive 



P(e) = Pi 
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(164) 



Given some initial values pi, et this relation describes, in the weak-field limit, the trajectory of the orbit on the (p, e) 
plane. Such curves for astrophysically relevant initial par ame ter s, ar e shown in Fig. (these curves remain essentially 



unchanged when the spin terms are retained in eqns. ( 138 ), ( 139 ) ). One feature that is immediately seen in the 



Newtonian-order inspiral is the absence o f the critical curve e = and the subsequent e > behaviour. This should 



not come as a a surprise, as expression (164) is not formally valid unless p 3> M. It is not safe to use weak-field 
approximations in strong field regimes. 



A simple way to make better predictions is by using the exact expressions (122), ( |123[ ) for p, e, but still employing 
the weak- field formulae ( |129| ), ( |130[ ) for the fluxes. The outcome of this trick is also shown in Fig. ^ for a set of 
orbits with initial parameters rp = 5,10,20M and ra = lO^M (this translates to e = 0.99999,0.99998,0.99996 and 
p = 10, 20, 40M respectively) for a = 0, 0.5M. We also considered a set of retrograde orbits with initial rp — 7 , 10, 20Af 



and the same apastron as before, and a — 0.99M spin. Note that these curves, like the ones given by eqn. (164), are 
shape-invariant with respect to the mass ratio as long as /i/M <C 1. It is rewarding to see that these new trajectories 
do show the existence of the e > region, and additionally are in good qualitive and quantitive agreement with the 
accurate numerical results |57|] . This is true as long as we do not attempt to evolve prograde orbits around rapidly 
spinning black holes, because the agreement quickly degrades when p becomes small. Essentially this approach takes 
into account the correct form of the potential, which is the main cause behind the change in sign of e, but the fact 
that the PN fluxes are increasingly inaccurate in strong-field regions precludes precise numerical agreement with the 
real trajectories. 

The new curves clearly predict a higher residual eccentricity as compared to the pure Newtonian curves. This 
approximate result strongly suggests that many astrophysically relevant inspiralling orbits will have a significant 
amount of eccentricity left when they are close to plunging and will therefore be likely to exhibit zoom-whirl behaviour. 
Our full numerical results cannot at present be used to reproduce complete trajectories, but Fig. ^, which displays 
arrows which are tangential to these trajectories at individual points, certainly shows that if significant eccentricity 
remains at p ^ 5M, that this eccentricity will not disappear in the last part of the inspiral before plunge. 

The work of Freitag |jl^ suggests that the initial periastra of scattered compact bodies which will eventually plunge 
into the black hole due to radiation reaction will be generally less than 40A/. However below that point the distribution 
of their periastra will be fairly flat, so that small initial periastra will be just as likely as large ones. As our figure 
shows, one expects, in the case of a Schwarzschild black hole, that bodies with initial rp > 20M will have e < 0.1 by 
the time of plunge. But for initial rp < 20M the final eccentricity will be e > 0.1 and can easily be as great as e ~ 0.7 
(see Fig. ^ or higher. For instance, for rp — lOM, the final eccentricity will be e '--^ 0.3 (see Fig. ^). We know that 
retrograde orbits will have less time to circularize and a longer "de-circularizing" time, so eccentricities in this case 
should be greater. This is clearly seen in Fig. ||. In the case of prograde orbits we should generally expect smaller 
eccentricities before plunging (as compared to orbits around Schwarzschild black holes) but still at a significant level. 
Again Fig. || suggests that the change in eccentricity will not be great, despite the longer circularizing and shorter 
de-circularizing times. Moreover, near extreme holes allow a wider range of initial periastra (for prograde orbits), 
even down to rp 2M and in these cases the residual eccentricity will be quite large. One concludes that eccentricity 
will play an important role in signal analysis for LISA. 

The technique outlined here has been recently applied in constructing approximate, radiation reaction-driven, 
inspirals of test-bodies in generic Kerr orbits (for which case only weak- field results are currently available [^). For 
more details see M. 
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FIG. 6. The radiative inspiral of a set of equatorial eccentric orbits with initial parameters (solid curves from left to 
right) rp = 5, 10, 20M and — lO^A'I, and for black hole spins a = 0, 0.5M. An additional set of retrograde orbits with 
Tp = 14, 20, 40M and same apastron, and for a — 0.99M is also shown (solid curves on the right side). The three dashed lines 
represent the separatrices for corresponding spins. Dotted curves are the Newtonian-order predictions, while the solid curves 
are the result of a more accurate calculation discussed in the main text. Note the significant qualitive difference between the 
two calculations at the vicinity of each separatrix. 



Motivated by the results discussed above, we would like to obtain simple estimates of the total amount of eccentricity 
gain and the number of orbits the particle will spend in the e > region. Such estimates may provide a useful guideline 
for assessing the observational importance of this phase. 

The numerical radiation reaction evolution-arrows in Fig. || show that the gradient de/dp grows (in the negative 
sense) monotonically as soon as the critical curve is cr ossed. Hence, the maximum (negative) value is attained exactly 



at the separatrix. The approximate expression (151) for de/dp is expected to be very accurate there. We can use 
this fixed gradient and extrapolate out to the critical curve, at some eccentricity e^, for a given eccentricity e/ at 
the separatrix. Then we have that 5e — Cf ~ Ci ~ {ps{ef) — Pc{ef)){de/ dp)p^(^^y This number should set an upper 
limit of the total increase in eccentricity. Results for some representative cases are given in Table From these 
numbers we deduce that, at best, there is a fractional increase of 5-50 % in eccentricity, the most favourable case 
being low-eccentricity retrograde orbits around rapidly spinning black holes. This gain decreases as we move upwards 
to larger final eccentricities (basically due to the shrinkage of the e > region). We therefore conclude that we should 
not expect any dramatic increase in eccentricity when the orbit is about to become unstable. 

A crude estimate on the number of orbits can be made by integrating ( |149| ) to find the time required to cross the 
e > region, 

t^^i r ^'^p (165) 

where we have factored out the angular momentum fiux as it can be taken as constant within the integration interval 
(and recovered by our numerical data). Similarly we have assumed a fixed eccentricity. The number of orbits is then 
calculated by dividing tc with a typical period Tr (or T^). We give some representative results in Ta ble |V^ . We need 
to emphasise that these numbers should be viewed only as order-of-magnitude estimates, as eqn. ( p_65| ) is a rough 
approximation. Nevertheless, we can still draw some reliable conclusions. For a small eccentricity our numbers are in 
agreement with existing results for nearly circular orbits around a Kerr black hole For e = 0.1 we should typically 
have a few thousands revolutions in the e > region around a a = 0.99M hole and for a mass ratio ji/M ^ 10~^. For 
the same parameters, but for retrograde orbits, there is an order of magnitude increase in the number of orbits. On 
the other hand, for all cases, there is an order of magnitude (or more) decrease as we move to eccentricity e = 0.5. 
Note that it is possible to have a small number of full orbits, but yet a significant number of azimuthal revolutions or 



"whirls" (see for example the a — 0.99, e — 0.5 case in Table VI) 
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C. Waveforms and fluxes from zoom- whirl orbits 



Let us now focus on the class of orbits w e hav e named zoom-whirl. As we have shown, orbits located near the 
separatrix should radiate in accordance with (147), as though they were nearly circular. In Table IV we list numerical 
fluxes for zoom-whirl or bits o f various eccentricities. It is clear from these results that as the separatrix is appr oached, 
the analytic prediction (147) is indeed confirmed. However, one has to be very cautious when applying (147) to the 
study of a real astrophysical, extreme mass ratio, binary system. As our data reveal, in the region where this relation 
is fractionally accurate at the level of ~ 10~^, the adiabaticity constraint (118) on the mass ratio is quite severe, 
typically n/M < lO^^ - iQ-s. 

The zoom- whirl orbits are of interest for future detection efforts because of the characteristic waveform they generate. 
In Figs. to ^ we show such waveforms (in particular we plot the quantity (/i/r)/i+ as a function of retarded time 
t — r*) for a couple of strong-field zoom-whirl orbits ( p = 2.11M and e — 0.7 , p — 1.7M and e = 0.3) and for a 
rapidly spinning black hole with a = 0.99M. Both orbits would evolve adiabatically if we consider a typical mass 
ratio ^/M ~ 10^^. The corresponding gravitational wave flux data can be found in Table VII. 

First we shall discuss the waveform as seen by an observer located on the black hole's equatorial plane, see Fig. |^ 
and top panel of Fig. ^ Clearly, these waveforms have a very distinct appearance. A rapidly oscillating, high 
amplitude signal is radiated during the whirling of the particle near periastron. In between these bursts one observes 
low-amplitude signals produced during the particle's zoom in and out from apastron. This contrast in amplitudes is 
greatest for larger eccentricities (compare Fig. |^ and Fig. An interesting feature of the equatorial waveform is the 
prominent high frequency ripples superimposed on the waveforms, associated with the higher multipoles components 
(£ = 3 and higher. The illustrated waveform includes all multipoles up to ^ = 18) of the wave. It is noteworthy 
that the high frequency features are prominent in both the zoom and the whirl parts of the orbit, although if they 
are solely the result of beaming we might expect that they would be purely a whirl feature, as the motion is fastest 
near periastron. However preliminary results from a time domain code written by one of us (KG) suggest that the 
high frequency features may be associated with quasi-normal modes of the black hole, which are excited by the high 
frequency emissions from the orbit. The quasi- normal mode ringing results in a continuous and time-delayed emission 
at these frequencies. 

We also show waveforms seen by an observer on the polar axis of the black hole, in the top panel of Fig. ^ and 
the bottom panel of Fig. ^ In this case both "plus" and "cross" polarisations are present (only h+ is nonzero for 
an equatorial observer) but we illustrate only /i+ because the "cross" waveform is the same except for a phase lag. 
The polar waveforms have the characteristic features of a high-amplitude, multi-cycle whirl part and a low-amplitude 
two cycle zoom part, but the high frequency features are absent. This suggests that the high frequency features are 
associated with beaming due to the rapid motion of the particle in the equatorial plane in the very strong field region, 
although the high-frequency features are not associated only with the strongest-field whirl part but are distributed 
throughout the whole cycle, including the weaker field (larger radius) zoom part. The polar waveforms are completely 
dominated by the quadrupole (£ = 2) emission. The £ = 3 and higher multipoles do not contribute significantly. In 
the equatorial waveforms the quadrupolar contribution is not much greater than that of the £ = 3 multipole and the 
fall off, in terms of the amplitude of hj^, for each subsequent multipolar waveform is slow. The transition between 
the polar and equatorial waveforms can by understood by looking at the waveform depicted in the bottom panel of 
Fig. H which corresponds to observation at angle 9 = 7r/4. 

One can get an idea of what is going on by looking at the ^-dependence of the energy flux from the system. In 
Fig. |l2| one sees that the m = 2 flux (dominated by the I = ra — 2 contribution) is concentrated somewhat towards the 
pole [9 — 0), but there is a strong shift towards the equator (0 = 7r/2) in the m = 3 flux, and further concentrations in 
that direction for each successively higher value of m. Therefore one sees the sort of beaming of the higher multipoles 
observed in the waveforms, although because of the dominance of the quadrupole the amplitude of the polar and 
equatorial waveforms is similar in this case. 

Next, in Fig. |lO|, we show the waveform from the retrograde zoom-whirl orbit p — lO.SAf, e — 0.5, retaining 
a — 0.99Af for the black hole spin. The familiar zoom-whirling pattern is clear also in this case. However, we do not 
see any prominent high-frequency structure in these waveforms as the contribution coming from higher multipoles 
is small (because the orbit does not reside in a very strong field regime). We are not surprised, in this case, that 
the waveform seen when observing from a point along the hole's polar axis (see bottom panel of Fig. |l^) is not very 
different. 

Although these zoom-whirl waveforms appear rather complex, and could present a problem for matched filtering 
data analysis techniques, one can take comfort in the fact that the number of harmonics of the spectrum which 
contribute significantly to the waveform and overall flux is not that great. Scott Hughes |Q has suggested that 
analysing individual "voices," which are monochromatic, of a complex wave may be an effective way to proceed in 
data analysis. Take the case of the zoom-whirl orbit with p = 2.11M, e = 0.7 and a = 0.99M. If we look at Fig. |lj. 
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which shows the entire spectrum (up to £ = 18) of this signal, we see that the number of individual harmonics, or 
voices, each corresponding to an instance of the three numbers £,m and k, which stand out are on the order of a dozen. 
Each of these voices has a rather simple waveform, as can be seen from Fig. ^ It is only their superposition which 
is complicated. Therefore following the chirp of individual voices as the orbit evolves may not be such a formidable 
computational task. 
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FIG. 7. The waveform produced by a particle in a zoom-whirl orbit with parameters p — 2.11M, e = 0.7. Specifically, we 

graph the quantity {fi/r)h+ (where r is the distance to the observation point, which is taken to be on the hole's equatorial 

plane) versus the retarded time t — r,(r) (in units of Af). We have set the black hole spin at a = 0.99M and included up 

to ^ = 18 multipoles in order to generate this figure. The orbital period is Tr — 236. 8M and A'^,. — 10.5. Note the very 

characteristic shape of the waveform, which is a periodic succession of high-amplitude/high- frequency parts (coming from the 

whirling motion of the particle near the periastron) and intervening low-amplitude/low frequency parts (from the zooming in 

and out motion). On the bottom panel, the same waveform is graphed over a shorter time interval, offering a clearer view of 

its rich structure. 
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FIG. 8. The same waveform as in Fig. [7| as seen by an observer along the hole's polar axis 6 — (top panel) and along 
6 = n/A (bottom panel). A comparison with the equatorial wave of Fig. Q reveals a substantial suppression of the high frequency 
features. This is a result of the fact that the wave's higher multipole components (which are responsible for the small-scale 
structure) are mainly "beamed" to directions close to the equatorial plane. 
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FIG. 9. The waveform generated by a particle in a zoom-whirl orbit with parameters p = 1.7M,e — 0.3 (we have again 
assumed a black hole spin a = 0.99M and £max ~ 17). The orbital period is Tr — 221.36A4" and the number of revolutions in 
one period is Nr — 12.3. The top and middle graphs show the signal seen by an equatorial observer, while the bottom graph 
corresponds to a polar observer. The same qualitive behaviour discussed in the caption of Fig. M is also evident here. 



33 



0.6 




-0.6 I ' ' ' 1 

500 1000 1500 2000 




500 1000 1500 2000 



FIG. 10. The waveforms generated by the retrograde zoom-whirl orbit p = 10. 5M, e = 0.5, a = 0.99Af (and Tr = 709M, 
Nr = 3.2) as viewed by an equatorial observer (top panel) and by a polar observer (bottom panel). Note the absence of any 
small-scale structure and the similar appearance of the wave from different viewing angles, a clear evidence of small contribution 
from high £ multipoles. 
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FIG. 11. The waveform associated with the individual harmonics (or "voices") £ = m = 2 and k — kmax = 10 (top), fc = 9 
(middle). The combination of the fc = 9 — 11 voices gives the waveform at the bottom which already shows some zoom- whirl 
behaviour. The total (. = m = 2 signal finally resembles the waveform shown at the top panel of Fig. 
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FIG. 12. This figure shows the angular dependence ol the energy flux for the orbit with p = 2.11M, e = 0.7 and a black 
hole with a = 0.99Af (see Figs. |^ and g). The x-axis shows the coordinate 9 in radians and the y-axis shows the rate of energy 
emitted into a angle 0.01 radians wide, as a fraction of the total energy emitted (at that multipole). Reading the different 
curves as they are peaked from left to right (where the left hand side of the graph corresponds to the pole of the black hole, 
and the right hand side to its equator) we have the antennae pattern for m=2, m=3, m=4, m=5, m=6, m=7, m=8 and at the 
extreme right m=18 (including the negative m contributions in each case). One notes that above m=2 the polar emission is 
greatly suppressed and the peak direction tends ever more towards the equator. In bold, one sees the curve for all multipoles 
at once (m= 2 to 18), but with the bins 0.1 radians wide, which peaks at around 9 = 7r/3. 



D. Spectra 



The final part of our numerical results concerns the harmonic decomposition of the gravitational radiation fluxes. 
Specifically, we examine the k distribution of the energy flux, at infinity and at the horizon, for a given multipole 
channel €,m (the angular momentum flux spectrum exhibits a similar behaviour). In all the figures in this Section, 
expect Fig. we plot Eimk versus k. 

To begin with, in Fig. |l^ we present the energy flux spectrum, at infinity, of the I = m = 2, 3, 4 multipoles for 
an orbit with parameters p = 2.11M, e = 0.7, and for spin a = 0.99Af. As we have already mentioned, for such a 
strong-field orbit the higher multipoles give significant contributions to the total flux. The spectrum itself is composed 
of a series of "humps" which grow in height as k increases, up to the point where the maximum harmonic is reached 
at fc = kmax, after which the spectrum rapidly fades away. This behaviour closely resembles the one found in the 
Schwarzschild case pl|], although in the present case the "humps" show a somewhat less regular structure. As we 
are not dealing with a very high eccentricity orbit, we flnd that k^ax is not very large. Moreover, the spectrum peak 
shifts to higher k values as ^ increases. 



In Fig. 15 we graph the corresponding horizon fluxes for the same orbit and same multipoles as in Fig. 13, Both 
infinity and horizon fluxes peak at a common kmax value. Since the selected orbit resides in a strong field regime, 
the horizon flux is a respectable fraction 10%) of the flux at infinity, and moreover, is predominantly negative i.e. 
superradiant, as we should expect from the discussion in the previous Section. Note that the relative contribution of 
the quadrupole channel is much larger in the case of horizon fluxes as compared to the infinity fluxes. 

The full energy flux spectrum, up to = 18, for the p = 2.11, e = 0.7 orbit is shown in Fig. |lj in terms of frequency 
rather that k. 

Spectra belonging to retrograde orbits appear similar to the prograde spectra with the difference that they peak at 
negative values of k. For a typical situation see Fig. |l^ for the retrograde orbit oi p = 10. 4M, e = 0.5 and black hole 
spin a = 0.99M. 

Before closing this Section, we should point out that only the harmonics around kmax (typically 10-30 of them) give 
a significant contribution to the total flux, for given m. This statement applies for all eccentricities we examined, 
and could have an important impact on the computation of waveforms and fluxes from highly eccentric orbits in which 
case kmax attains very large values. For example, if one could find by other means (for example a time-domain code, 
see discussion in the next Section) the location of the main peak, then the calculation of the surrounding harmonics 
should give an sufficiently accurate result for the flux. This strategy is far more economic, from a computational point 
of view, than calculating all the harmonics between fc = and kmax- 
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FIG. 13. The distribution of the energy flux (in units of fX^ /M^) in terms of the k harmonic number for an orbit with 
p = 2.11M, e = 0.7 around an o = 0.99M Kerr black hole. Prom top to bottom we display the £ = m = 2, 3, 4, 5 
multipoles. Note the significant contribution of the higher multipole channels to the total flux. The relevant orbital frequencies 
are Mfir = 0.02653 and MOw, = 0.2791. 
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FIG. 14. This figure siiows all of the main peaks of the spectrum emitted by the orbit p — 2.11M, e = 0.7 for a black hole 
with a = 0.99M (see Figs. |^ and ^ for the waveform associated with this orbit). The flux of total energy emitted per unti 
time, on the y-axis, and the frequency, on the x-axis, are both given in the geometrised units of this paper (in which 5/iS is 
approximately unity, if the black hole has mass of lO^M©). The main peaks in each multipole are easily read from left to right as 
(£ = m^2,k = 10); (^ = m = 3, fc = 15); {£ = m^4;k = 21); (^ = m = 5, fc = 26); {£ = m = 6,k = 31); (^ = m = 7, fc = 37) 
{£ = m = 8,k = 42);{£ = m = 9,k = 47); (£ = m = 10, fc = 53); {£ = m = ll,k = 58); (^ = m = 12, fc = 63) 
(^ = m = 13, fc = 69); (^ = m = 14, fc = 74); (^ = m = 15, fc = 80); (£ = m = 16, fc = 84); (f = m = 17, fc = 90) 
{£ = m = 18,fc = 94). 
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FIG. 15. The horizon energy flux fc-spectrum for the orbit p = 2.11M, e = 0.7, a = 0.99M. As in the previous figure, we 
graph, from top to bottom, the multipoles £ = m = 2, 3, 4. Note that both infinity and horizon fiuxes peak at the same k 
harmonic, for given £, m. The latter spectrum, however, is strongly dominated by the quadrupole channel (which is roughly 
10% of the fiux at infinity) while the higher multipoles quickly fade away. The fact that the horizon flux takes, almost entirely, 
negative values means that it represents superradiant radiation. 
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FIG. 16. The £ = m — 2 (top) and £ — m — 3 (bottom) spectra for the p = 10.4, e = 0.5 retrograde orbit. The black 
hole spin is o = 0.99M . Note the location of the maximum at some negative A; value and the dominance of the quadrupole 
component over the octapole (and higher) component. 



VI. CONCLUDING DISCUSSION 



Wc have examined gravitational waves from and radiation reaction of equatorial orbits of particles in the last 
stages of inspiral around a central black hole. We expect that such orbits will often have moderate eccentricities, 
not circular but significantly less than e = 1. As these orbits approach the point at which they will plunge into 
the black hole two things will happen. The first is that the orbital eccentricity will begin to increase, having been 
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decreasing throughout the preceding inspiral. For retrograde orbits this eccentricity increasing phase will last for 
many cycles, while for prograde orbits, especially ones which retain fairly high eccentricities, this phase will be fairly 
brief (in the case of rapidly rotating black holes). We have found a ~ 10% fractional increase in eccentricity, for 
the most favourable situations. Secondly, as the orbit draws closer to the unstable region it will tend to spend a 
longer and longer portion of each orbital period near periastron. For reasonably eccentric orbits this will promote 
the "zoom-whirl" behaviour we have described, with its characteristic waveform. Although extreme cases involving 
high eccentricities will presumably be rare, it is still likely that LISA will be seeing signals with 5—10 whirls, and 
therefore 10 — 20 gravitational wave cycles between each apastron. Therefore templates of such signals will be very 
important. It is worth noting that because the orbital radius changes very little during these whirls, the whirl part of 
the waveform looks, in many respects, like a near circular waveform with harmonics of a single dominant frequency. 
This means, for instance, that for orbits with many whirls, the ratio of the fluxes of energy and angular momentum 
is close to (recall the celebrated relation for circular orbits E/L = fi^). 

Considering inclined orbits (not confined to the equatorial plane) for a moment, we expect that zoom- whirl orbits 
will be found in cases with small inclination angles, because higher inclination angles feature greater plunge radii and 
the particle cannot ever have a very small periastron radius (as we see when the orbit is retrograde, which corresponds 
to the largest possible inclination angle). When the test-particle is in an inclined and eccentric orbit, close to crossing 
some separatrix of bound stable orbits, it will spend a considerable amount of time moving in a quasi-circular non- 
equatorial trajectory close to the periastron. It is plausible, by extrapolating the results of the present work, that 
during the "whirling" stage, the orbital energy, angular momentum and Carter constant will approximately evolve in 
such a way as if the orbit was circular and inclined. Assuming that most of the radiated energy, angular momentum 
and Carter constant are generated near the periastron (which in the present case is located in a strong field regime) 
then one might be able to estimate the, otherwise elusive, rate of change of the Carter constant for a near-to-plunge 
generic Kerr orbit. Such information could provide a very useful test for the recently adopted assumption of obtaining 
the rate of change of the Carter constant by keeping the orbital inclination angle fixed during the inspiral [ p7| . 

Waveforms from prograde orbits residing close to the horizon of the black hole will feature significant high frequency 
components when seen from a position on the equatorial plane of the system (i.e. when the orbit is observed "edge 
on"). This seems to be due to beaming, resulting from the rapid motion of the orbiting particle along the line of sight 
of the observer. When observed from on or near the polar axis the waveform is largely quadrupolar, dominated by a 
single nearly circular frequency. A glance at the equatorial zoom-whirl waveforms presented here suggests that they 
can be very complex and not necessarily amenable to matched filtering methods following the full wavetrain. But 
recall that a successful source identification may already have been made during the previous year of observation by 
LISA and from the source parameters deduced during this period it may be possible to search for the whirl parts of 
the waveform individually. On the other hand, it is worth noting that this high frequency structure could make it 
possible for LISA to detect late inspiral signals from very large black holes, with masses above IO^Mq, which would 
otherwise be too low frequency (in the low multipole parts of the waveform) for detection. By contrast, if we are 
looking down on the system from the pole then the signal is much "cleaner" , without much contribution beyond the 
i = m = 2 multipole. Obviously non-equatorial motion will introduce further harmonics and, as suggested in |3^ ] 
it may prove more useful to examine different harmonics or "voices" of the signal separately, rather than trying to 
model the entire complex signal as one template. 

It has been recently realised that orbits during the long inspiral phase, before the radial orbital frequency is 
within the LISA waveband, will, in principle, emit detectable radiation. Recall that the periastron is rather close, 

< 20M , so that the frequency of the cycles in the whirl part of the waveform does fall within the LISA waveband. 
As the time between successive bursts (equal to the radial orbital period) will be very long (typically up to a century 
or even more) these bursts would not be expected to be detectable in practice. However, as the number of objects 
in this long inspiral stage will be rather great, one does have to take them into account as a background noise which 
will tend to hinder LISA's efforts to detect signals from those objects in the last stages of their inspiral. In fact, one 
could in principle expect to have proper zoom-whirl orbits even at this stage, provided the central hole is rapidly 
spinning and the orbit is retrograde. Then, as we have shown, the corresponding separatrix translates to a minimum 
periastron value, nsbo ~ 6Af , which falls within the expected periastron distribution. 

In a follow-up paper ||33|] we will study orbits with larger radii and larger eccentricities. In particular we are 
planning to produce (using the method of which was applied for inspiralling circular inclined orbits) the full 
inspiral trajectory and the resulting waveform for bound equatorial orbits that could be of importance for LISA. It 
is of great interest to (a) produce full inspiral waveforms for the kind of orbits described above and (b) give a good 
estimate of the total inspiral time, and the residual eccentricity (taking under consideration effects like the sign reversal 
of e) just before plunging. However, the calculation of fluxes and waveforms produced by bodies in e « 1 orbits, is 
currently ranked as a difficult task. When pursued in the frequency domain (as it was the case in the present work) 
one has to calculate an enormous number of individual harmonics, as it is obvious from the spectra we presented. 
Moreover, the numerical computation of the integrals (Isl) is poorly convergent, a manifestation of the fact that the 
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source term in the Teukolsky equation ( p4| ) diverges as r — > oo (that is, when the orbit tends to become parabohc 
e — *■ 1). One way to cure this pathology would be to work with the inhomogeneous Sasaki- Nakamura equation ]3^ ] 
which has a well behaved source term (in the sense that it decays at spatial infinity). However, the first difficulty 
outlined above will still be present. A possible way to overcome it could be the calculation of the waveform/fluxes 
directly in the time domain by evolving the time-dependent Teukolsky equation without resorting to any separation 
of variables apart from (p. Conceivably, the required numerical code could be based on the Teukolsky codes used to 
study the dynamics of scalar and gravitational perturbations in a Kerr background metric [ p8[ , see | p9| for a report 
on such an attempt. Looking further ahead, such time-domain codes could be the only practical tool for computing 
the waveform/fluxes generated by bodies orbiting non-black hole massive compact objects (in which case there is no 
known Teukolsky- like separable wave equation). We are currently working along both of these directions. 
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VIII. TABLES 
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a = 0.5A/ 


a = 0.99Af 


a = -0.99A'/ 


0.10 


4.377 (4.71) 


1.516 (1.59) 


9.266 (10.03) 


0.20 


4.526 (4.77) 


1.595 (1.64) 


9.552 (10.12) 


0.30 


4.679 (4.85) 


1.685 (1.71) 


9.830 (10.24) 


0.40 


4.836 (4.96) 


1.782 (1.79) 


10.102 (10.40) 


0.50 


4.996 (5.08) 


1.883 (1.89) 


10.367 (10.58) 


0.60 


5.158 


1.988 


10.627 


0.70 


5.323 


2.094 


10.882 


0.80 


5.490 


2.201 


11.133 


0.90 


5.658 


2.310 


11.380 


LOO 


5.828 


2.420 


11.623 



TABLE I. The separatrix ps and the critical value Pcrit where e = (in 
parentheses, accurate to the decimals shown) for a variety of eccentrici- 
ties and for three different black hole spins, a — 0.5AI, a = 0.99Af and 
a = — 0.99Af (retrograde orbits). 
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<? 


rcrtt/M 


corrected value 


-0.9 


9.64 


9.74 


-0.5 


8.37 


8.43 


u.u 


D.Do 


D.Do 


0.5 


4.70 


4.69 


0.7 


3.76 


3.75 


0.9 


2.56 


2.54 


0.95 


2.03 


2.11 


.99 


1.47 


1.55 


1.0 


1.0 


1.0 



TABLE II. The position of the critical radius, rcrit 
in units of M, for different black hole spins a and zero 
eccentricity. The parameter q — a/M is defined here 
to be negative for retrograde orbits and positive for 
prograde orbits. This table is provided as an erratum 
to Table I of Ref. which was incorrect due to a 
bug in the part of the code calculating the fluxes of 
energy and angular momentum radiated to the black 
hole horizon. This data was produced using the cor- 
rected code from the previous paper, rather than with 
the code of the present paper. 



a/M 


e 


p/M 


(M//i)'ijgW- 


0.95 





10.015 


4.966452E-05 
(4.966247E-05) 


0.95 





40.795 


5.277469E-08 
(5.277415E-08) 


0.95 





200.698 


1.933592E-11 
(1.933573E-11) 


0.00 


0.7641 


8.754 


1.57132E-04 
(1.57131E-04) 


0.00 


0.7446 


13.198 


1.43629E-05 
(1.43632E-05) 


0.90 


0.3731 


12.152 


2.3570E-05 
(2.3893E-05) 


0.90 


0.5634 


50.513 


2.1211E-08 
(2.1192E-08) 


0.30 


0.6519 


19.969 


2.1654E-06 
(2.1375E-06) 



TABLE III. Comparing results from our radiation reaction code with 
existing results found in the literature [42], |5^], jsi] (data in parenthe- 
ses). We find excellent agreement (at the predicted level) for equatorial 
circular Kerr and eccentric Schwarzschild orbits. On the other hand, 
there is seems to be a ~ 1% disagreement with Shibata's results for 
equatorial eccentric orbits. 
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a/M 


e 


p/M 


(-^gW/-^'gW)/^</> 


0.50 


0.3 


4.70 


1.071 


0.50 


0.4 


4.90 


1.128 


0.99 


0.3 


1.70 


0.984 


0.99 


0.3 


1.80 


0.896 


0.99 


0.4 


1.80 


0.976 


0.99 


0.7 


2.11 


1.138 



TABLE IV. Examining the validity of the approximate, 
near-separatrix, formula E = Q.^L for various zoom-whirl orbits and 
for two black hole spins. Here, only the fluxes at infinity have been con- 
sidered (the horizon fluxes yield similar results). Typically, this relation 
is found to be accurate to fractional accuracy 10-^-10-2. In such cases 
the orbit is so close to the separatrix as to require fi/M <C 10~^ — 10~^ 
for adiabaticity to hold. 



a/M 


e/ 




M{de/dp)p^ 


Se/ci 


0.50 


0.1 


0.086 


-0.0414 


0.16 


0.50 


0.3 


0.28 


-0.1342 


0.073 


0.99 


0.1 


0.086 


-0.1555 


0.16 


0.99 


0.3 


0.29 


-0.2186 


0.019 


-0.99 


0.1 


0.066 


-0.0439 


0.52 


-0.99 


0.3 


0.28 


-0.0546 


0.083 



TABLE V. Upper limits on the total eccentricity gain close to the 
separatrix, for given flnal values e/ for the eccentricity. Using the gradi- 
ent de/dp at the separatrix we extrapolate to the critical curve. In this 
way we obtain the eccentricity d. 



a/M 


e 


{^i/M''% 


Tr/M 


n/M 


(p/M)Nr 


{li/M)N4, 


0.99 


0.1 


0.051 


216.48 


18.03 


2.3E-4 


2.8E-3 


0.99 


0.5 


0.0018 


276.84 


18.02 


6.7E-6 


l.lE-4 


-0.99 


0.1 


5.6 


651.48 


181.02 


8.6E-3 


3.1E-2 


-0.99 


0.5 


0.79 


718.32 


218.14 


l.lE-3 


3.6E-3 



TABLE VI. Approximate data for the number of orbits in the 
e > regime. The required crossing time is tc and we have defined 
Nr,4i = tc/Tr,^. We have calculated the periods Tr,^. aip = {ps+Pcrit)/^. 



a 

e 


= 0.5M 

p/M 


(M/fifE^^. 


{M/iifEBw 




{M/t?)L^^ 


(M/m)p 




0.10 


4.60 


2.88029E-3 


-6.41673E-6 


2.88686E-2 


-6.47058E-5 


-5.79612E-1 


9.61181E-3 


0.10 


5.00 


1.81736E-3 


-4.03021E-6 


2.06724E-2 


-4.58643E-5 


-2.37079E-1 


-3.73062E-3 


0.10 


6.00 


7.10665E-4 


-1.27547E-6 


1.05539E-2 


-1.88243E-5 


-8.45033E-2 


-2.00595E-3 


0.20 


4.70 


3.11812E-3 


-5.71886E-6 


2.96692E-2 


-5.63759E-5 


-5.65181E-1 


1.04172E-2 


0.20 


5.00 


2.09142E-3 


-4.27185E-6 


2.21228E-2 


-4.58884E-5 


-2.56971E-1 


-7.19373E-3 


0.20 


6.00 


7.78541E-4 


-1.43776E-6 


1.08487E-2 


-1.97306E-5 


-8.53506E-2 


-3.97940E-3 


0.30 


4.70 


4.96241E-3 


-2.34257E-6 


4.07599E-2 


-2.62670E-5 


-2.55152E-0 


2.64021E-1 


0.30 


5.00 


2.60439E-3 


-3.79849E-6 


2.48384E-2 


-3.98856E-5 


-3.03745E-1 


-9.63732E-3 


0.30 


6.00 


8.88282E-4 


-1.63235E-6 


1.12981E-2 


-2.06383E-5 


-8.65523E-2 


-5.87677E-3 


0.40 


4.90 


4.52598E-3 


3.00302E-6 


3.62936E-2 


9.69843E-6 


-9.12063E-1 


3.28050E-2 


0.40 


5.00 


3.53043E-3 


-5.18531E-9 


2.98097E-2 


-1.10071E-5 


-4.41342E-1 


-6.78246E-3 


0.40 


6.00 


1.03261E-3 


-1.68472E-6 


1.18257E-2 


-2.03399E-5 


-8.77285E-2 


-7.62622E-3 
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0.50 


5.10 


4.21594E-3 


9.19757E-6 


3.26383E-2 


5.01972E-5 


-5.47629E-1 


-4.96045E-3 


0.50 


5.50 


2.11797E-3 


6.77254E-8 


1.89546E-2 


-9.89420E-6 


-1.60976E-1 


-1.41868E-2 


0.50 


6.00 


1.19638E-3 


-1.22374E-6 


1.22973E-2 


-1.65642E-5 


-8.82324E-2 


-9.10645E-3 




a - 


= 0.99M 














e 


p/M 


(M/M)^£gW 




(M//i^)iS=^. 


(M/M^)igw- 


{M/^^)p 


(MV/i)e 


0.10 


1.55 


9.26325E-2 


-7.85155E-3 


2.63428E-1 


-2.23134E-2 


-1.51950E-0 


1.74361E-1 


0.10 


2.00 


4.72325E-2 


-3.16550E-3 


1.77532E-1 


-1.18650E-2 


-4.73445E-1 


-3.72486E-2 


0.10 


3.00 


1.12400E-2 


-4.14404E-4 


6.8334 7E-2 


-2.50408E-3 


-2.26534E-1 


-1.26041E-2 


0.20 


1.62 


9.30011E-2 


-7.62204E-3 


2.60464E-1 


-2.13080E-2 


-1.46506E-0 


1.57196E-1 


0.20 


2.00 


5.06654E-2 


-3.43868E-3 


1.82214E-1 


-1.22541E-2 


-4.74856E-1 


-7.36986E-2 


0.20 


3.00 


1.19893E-2 


-4.59376E-4 


6.94427E-2 


-2.60733E-3 


-2.25311E-1 


-2.48251E-2 


0.30 


1.70 


9.56364E-2 


-7.52010E-3 


2.59798E-1 


-2.03993E-2 


-1.64242E-0 


1.74157E-1 


n on 
O.oO 


2.00 


C CO /I 1 OT^ O 

0.6o4izrj-z 


-6.oio26rj-6 


1 onn /1 1 1 


-i.20D4yrj-2 


A '7000 1 "CP 1 
-4. /OOoili-i 


1 AOKIOT^ 1 

-i.U8oizrj-i 


0.30 


3.00 


1.31541E-2 


-5.29365E-4 


7.10070E-2 


-2.76123E-3 


-2.22638E-1 


-3.62186E-2 


0.40 


1.80 


9.53548E-2 


-7.06881E-3 


2.55921E-1 


-1.89132E-2 


-1.28174E-0 


-5.85767E-3 


0.40 


2.00 


6.42861E-2 


-4.43924E-3 


2.00778E-1 


-1.36388E-2 


-4.90607E-1 


-1.40712E-1 


0.40 


3.00 


1.45880E-2 


-6.16038E-4 


7.25434E-2 


-2.93630E-3 


-2.17455E-1 


-4.61972E-2 


0.50 


2.00 


7.50848E-2 


-5.10737E-3 


2.15885E-1 


-1.45086E-2 


-5.33508E-1 


-1.68200E-1 


0.50 


2.50 


3.29957E-2 


-1.83305E-3 


1.22205E-1 


-6.59015E-3 


-3.02830E-1 


-9.41048E-2 


0.50 


3.00 


1.60427E-2 


-7.05424E-4 


7.32601E-2 


-3.08563E-3 


-2.08134E-1 


-5.39259E-2 


0.70 


2.11 


9.29845E-2 


-5.01552E-3 


2.39102E-1 


-1.30762E-2 


-9.57156E-1 


-1.64115E-1 



a = -0.99M 



e 


p/M 




{M/^jifEBw 










0.10 


9.5 


1.22528E-4 


1.50991E-6 


-3.31424E- 


3 


-3.98335E- 


5 


-1.93846E- 


1 


4.94557E-3 


0.10 


11.0 


4.99506E-5 


3.39590E-7 


-1.72497E- 


3 


-1.13847E- 


5 


-3.05501E- 


2 


-2.63944E-4 


0.20 


9.7 


1.40484E-4 


2.21408E-6 


-3.53943E- 


3 


-5.25871E- 


5 


-2.20308E- 


1 


7.70198E-3 


0.20 


11.0 


5.63927E-5 


5.01279E-7 


-1.80663E- 


3 


-1.47518E- 


5 


-3.17662E- 


2 


-5.24602E-4 


0.30 


10.0 


1.49711E-4 


2.90957E-6 


-3.53885E- 


3 


-6.33565E- 


5 


-1.70530E- 


1 


4.06284E-3 


0.30 


11.0 


6.74024E-5 


8.29938E-7 


-1.94163E- 


3 


-2.11328E- 


5 


-3.40152E- 


2 


-7.76820E-4 


0.40 


10.3 


1.57135E-4 


3.72816E-6 


-3.47105E- 


3 


-7.49426E- 


5 


-1.33052E- 


1 


1.49356E-3 


0.40 


11.0 


8.33663E-5 


1.42912E-6 


-2.12770E- 


3 


-3.18078E- 


5 


-3.77928E- 


2 


-1.00954E-3 


0.50 


10.4 


2.46763E-4 


8.15567E-6 


-4.73640E- 


3 


-1.46459E- 


4 


-5.62681E- 


1 


1.91167E-2 


0.50 


11.0 


1.04907E-4 


2.47885E-6 


-2.36298E- 


3 


-4.88973E- 


5 


-4.44577E- 


2 


-1.19582E-3 



TABLE VII. Numerical data for the rate of change, under radiation reaction, of E, L (separately for infinity and the horizon) 
and p,e (total amount) for a selection of strong-field orbits and for two black hole spins, a = 0.5Af (top), a — 0.99M (middle) 
and a = — 0.99Af (bottom). Most of these data were used to generate the vectors in Fig. In the computations we have used 

^max = 10 - 17. 
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APPENDIX A: FUNCTIONS THAT APPEAR IN THE SOLUTION FOR . 

The quantity x — Lz — aE satisfies the quartic equation, 

e)a;4 + N{p, e)x^ + C{p, e) = (Al) 



where 



F{p, e) = - 2M(3 + e2)p2 + M^{3 + e^fp - ^Ma'il - e^f] (A2) 



N{p, e) = -[-A-V + (M^(3 + e^) - a^)^ - Ma^{l + Se^)] (A3) 
P 

C(p, e) = (a^ - Mpf (A4) 
Defining the discriminant 



A^(p, e) = _ 4^^- = - 4Afp3 + 2{2A/2(i _ + ^a^^ ^ ^2)1^2 „ ^Ma^{l - e^)p + a4(i _ e^)^] (A5) 



p 

the solution for is, 

2 -iVTAy' .... 
^ = (A6) 

where the upper (fower) sign corresponds to prograde (retrograde) motion. 

APPENDIX B: FUNDAMENTAL FREQUENCIES FOR BOUND EQUATORIAL ORBITS IN KERR 

GEOMETRY. 

The fact that the function r(t) is periodic in time (with period T-^) imphes that the function 

dd) aT + Ax 



(Bl) 



dt (r2 + a^)T + axA 
is also periodic and with the same period. Hence, it can be expanded in a Fourier series, 

^ = E /^ke-'^^^* (B2) 



dt 

k— — oo 



By integrating this relation we get. 



(l,{t) =(3ot + Y^ Cke-*'^^^'* + (const) (B3) 



where Ck = i(3k/kflr. This expression clearly shows that the (/)-motion is of "rotation" type [Q. Making use of the 
integration constant to define the term cq, we have 

</)(0-/3oi = Ecke-''''''* (B4) 

k 

That is, the function F{t) = (j>{t) — (i^t is periodic with period equal to T,.. As we have defined = 0(< + T,) — (j){t) 
we find that 

/3o - ^ ^ (B5) 
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APPENDIX C: POTENTIALS OF THE SASAKI-NAKAMURA EQUATION 



We give an explicit listing of the potentials F{r),U{r) that appear in the Sasaki-Nakamura equation (10£) 



F{t) = 



U{r) = + G{r f + ~ F{r)G{r) 



(r2 + a2) 



where the function ri(r) is given by 



with the following coefficients 



ri{r) = Co + ci/r + C2/r^ + Ca/r^ + C4/r* 



Co = -f 2icjM + A(A + 2) - 12aw(aw - to) 
ci = 8ia[3acj — A(aa; — m)] 

C2 = -24zaAf {au - to) + f 2a2[l - 2(aw - m)^] 
C3 = 2Ua^{atjj - m) - 24Ma^ 



C4 = f2a*. 

In addition, the functions G{r) and Ui{r) are 

2(r-M) rA 



f/i(r) = V{r) + 



A2 



2a 



1^ ,r\ V,r f P,r 



with 



i(3K 



ZiK,r + A + ^ 



6A 
T 



r 



(CI) 
(C2) 

(C3) 



(C4) 
(C5) 
(C6) 
(C7) 
(C8) 



(C9) 

(CIO) 
(CII) 

(CI2) 



(3 = 2A{-iK + M ~ 2A/r) 
The functions K{r), V{r) appear in the Teukolsky equation (pl|). 
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